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Abstract 


A non-equilibrium stage model has been developed for the steady state and dynamic 
simulation of a countercurrent multicomponent liquid-liquid extraction column. Dis- 
persed phase behavior described by means of drop population model, allows to split 
the drop spectrum into different drop classes interacting with bulk continuous phase. 
Mass balances are then written for each drop class and bulk continuous phase in a 
stage. Rigorous transport relations using interactive mass transfer coefficients have 
been used to explain interphase mass transfer in extraction column. This model to 
incorporates effects such as reverse and osmotic diffusion which remain unaccounted 
for by simple transport relations. Drop breakage and coalescence have been consid- 
ered with mass transfer in this study; generation term Uij k accounts for these in mass 
balance equations. 

Newton- Raphson (N-R) method has been used to solve the set of non-linear equa- 
tions for steady state. DASSL (Differential Algebraic System SoLver) with N-R 
method has been implemented to solve differential-algebraic equations (DAE) for 
dynamic simulation. DASSL has various advantages over other methods such as 
Runge-Kutte and Gear’s method apart from being simple and easy to understand. 
Computer program (simulation package) of about 12,000 lines of “C” code divided 
into 19 files has also been developed as a part of this thesis. Model developed in this 
thesis is general and applicable to all types of multicomponent counter-current extrac- 
tion column. In this work the model has been tested for unagitated sieve tray column. 
Two extraction system, Aromatic extraction by Sulfolane and Acetone extraction by 
Toluene have been taken for this study. Voluminous results generated by the sim- 
ulator have been thus analyzed this covers study of concentration profile, dynamic 
holdup, flow and velocity profiles under various feed conditions and perturbations. 
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Chapter 1 
Introduction 


1.1 General 

Chemical process industries, often need to separate the components of homogeneous 
liquid mixture. Out of several approaches, liquid-liquid extraction (LLX) commonly- 
known as solvent extraction, is a technique for separating liquid components of a so- 
lution by distribution between two liquid phases. With the recent emphasis on energy 
conservation in process industries, extraction today is a choice preferred over distil- 
lation especially for separation of dilute liquid mixture in which distillation requires 
enormous amount of energy. Mixtures of heat sensitive components can be separated 
at low temperature using suitable solvent by LLX. 

Solvent used in the extraction process should be immiscible or partially miscible 
with one of the component of mixture in order to facilitate the separation of liquid 
phase. Liquid-liquid extraction is now being adopted as a more economic alternative 
to other separation processes and has found immense application in the separation 
of: 

• Solution of components having low relative volatility, especially when vacuum 
distillation is expensive. 

• Solution of close boiling point and azeotrope forming components. 

• Separation of solute from mixture when evaporation is impractical. 

• Solutes of heat sensitive components such as antibiotics. 

Over a period of time organic, inorganic and nuclear industries developed interest in 
LLX operation. V arious processes such as aromatic extraction, acetic acid extraction 
from water, Uranium and other hydrometallurgical extraction developed commercially 
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have been in operation for a long time. In order to understand these processes thor- 
oughly and generalize results from one process to another, theoretical research is 
necessary. Modeling of separation processes has now become an integral part of 
chemical process development and scale up process. With the help of theoretical 
models, process behavior could be visualized a priori and important decisions can be 
arrived at before erecting actual equipment at site. This not only saves money but 
provides better insight of the process also. Even though much of the work with which 
chemical engineers are involved, concerns the steady state behavior of process (which 
are quite useful in equipment design decisions) the unsteady state characteristics of 
a process are also important in a certain class of problems such as : 

Startup problem in which it is desired to predict the rate at which a process pro- 
ceeds to steady state with given initial state. 

Control problem in which information is desired on the open loop response of 
the controlled variables in process due to uncontrolled disturbances in certain 
input stream. With the advent of fast computing machines, these models can be 
tailored for on line process simulation and control (such as feed forward control). 
Prior study could be carried out on process simulator to design control system 
for the process plant which is yet to be erected. 

Training simulator Training of operators for process, without having recourse to 
real equipments, could be facilitated by a simulation model, which imitates the 
real process. 

Present thesis is devoted to the development of both steady-state and dynamic 
simulation model for multicomponent countercurrent liquid-liquid extraction process. 


1.2 Literature review 

Modeling of extraction process 

One of the simplest way to model any chemical separation process such as extraction, 
distillation, absorption for a staged process is to assume equilibrium state at each 
stage. Such an approach is rather impractical, especially in extraction process where 
extraction efficiency may be as low as 20-25%, Laddha [1]. Hence an alternate, rate 
based approach is proposed to model the mass transfer behavior of an extraction 
process. Several attempts have been made in this direction such as by Ricker [2], 
Spencer [3] and others. Apart from non-equilibrium, staged processes also experience 
the phenomena of backmixing or axial dispersion. Baird [4], Ricker [2], Spencer [3] 
have dealt with this phenomena by assuming some fraction of main flow as the back- 
mixing flow. Recent success of rigorous models based on coupled fluxes in distillation 
(Krishnamurthy [5, 6]), encourages the use of such models for extraction process mod- 
eling. Zimmermann’s [7] proposed model for extraction allows rigorous calculation 
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of mass transfer coefficients incorporating interaction effect. Zimmermann extended 
the non-equilibrium approach with drop population balance and rigorous mass trans- 
fer coefficient calculation. Earlier work by Khani [8] using the population balance 
approach has been formulated as set of partial differential equations. However mass 
transfer aspects were not discussed in detail. 

Modeling of_ drop breakage arui coalescence 
Extraction involves two distinct phases - Continuous and Dispersed. Dispersed phase 
is described by the droplets formed in the process. During the past three decades 
a need has been felt to incorporate drop population approach. Review article by 
Ramakrishna [9] described the status of population balance vis a vis other aspects 
of modeling of liquid dispersion. In the direction of liquid-liquid dispersion one of 
the pioneering work was by Coulaloglou and Tavlarides [10]. Tavlarides proposed 
a phenomenological model to describe drop breakage and coalescence in a turbu- 
lently agitated liquid-liquid dispersion. In last 3 decades several workers have extended 
this approach to modeling of extraction process. Sovova [11] worked out the same 
model for batch stirred extraction process and discretized it based on the real stages. 
Sovova [12] used the same discrete model and applied drop population approach for a 
vibrating plate extractor Hsia and Tavlarides [13] used a similar approach for stirred 
tank. Hsia and Tavlarides [14] described drop breakage, coalescence and micromix- 
ing in stirred tank. Rubio [15] discusses about the various possible theoretical drop 
spectrum distribution in a Wirz extraction column and fits them to a real extraction 
operation. Recent work by Tsouris and Tavlarides [16, 17, 18] further discusses in 
great deal about the model development for drop breakage and coalescence based on 
earlier work of Coulaloglou [10]. The model has been analyzed for drop breakage and 
coalescence behavior and control of dynamic holdup along the multistage extraction 
column. 

Mass transfer and axial dispersion 

Above works on drop population behavior, incorporate only drop population model 
and say very little about mass transfer effect on extraction process. Zimmermann [7] 
synthesized both drop population and non-equilibrium approach and proposed a de- 
tailed mathematical model for extraction column at steady state. The approach relies 
a Poly dispersed phase model of Jiricny [19] according to which dispersed phase having 
a drop spectrum is assumed as several separate phases (plugs of dispersed phase - 
each plug contains drops of same volume), interacting with continuous bulk phase. 
Backmixing is assumed in form of axial dispersion between adjacent stages, which is 
represented by axial dispersion coefficients D c ,D d . Rigorous transport relations as 
used by Krishnamurthy [5, 6] for distillation and absorption are used for extraction 
process also. These transport relations have been independently solved by lineariz- 
ing relationships Stewart and Prober [20], utlizing the film theory of Taylor [21] and 
Smith [22]. In this thesis these relations are the part of a simulation model thus 
solved simultaneously with other equations. 

Equilibrium relations 
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Every extraction model uses equilibrium relations in form of yi = KiXi, where distri- 
bution ratio K{ is one of the tedious factors in modeling of extraction column. For 
organic systems especially used in refinery and petrochemical complexes, UNIFAC 
method is one of the most relied and easy to use method available in literature. The 
method based on group contribution has been proposed by Fredenslund and Praus- 
nitz [23] for non-ideal liquid mixtures. Pioneering work by Magnussen [24], reports 
UNIFAC LLE parameters for 32 main groups and 57 sub-groups in the temperature 
range of 10 — 40°C. LLE data for specific systems has been reported in DECHEMA , 
Chemistry Data Series part I, II, III [25]. 

A modified UNIFAC group contribution model Larsen [26] also allows to incor- 
porate temperature effect, but published data on this model is scarce. Several other 
papers published on specific systems such as LLE phase equilibrium for dearomati- 
zation of ATF (Aviation turbine fuel) fraction by Mehrotra, Garg and Chopra [27]. 
Mukhopadhyay [28] published UNIFAC parameters for aromatic extraction system 
for high temperature also. Hooper [29] published UNIFAC data in the temperature 
range of 20 — 250°C for some water organic liquid system. 

Empirical correlations 

Several correlations used for extraction columns such as characteristics velocity and 
binary mass transfer coefficient have been refered to from the design book by Lad- 
dha [1] and handbooks by Baird [4] and Rousseau [30]. 

Solution of algebraic and differential -algebraic equations system 
For steady state simulation Block-tridiagonal equations sets are solved by Newton- 
Raphson (N-R) method and Block-thomas algorithm. However for dynamic sim- 
ulation a set of Differential-Algebraic Equation (DAE) is formed. Review paper by 
Gani [31] and Gani [32, 33] applied the solution for DAE system. Several DAE system 
solver are proposed by authors such as by Gear’s [34, 35] and variety of Runge-Kutta 
methods which are popular for solving DAE. In recent years DAE solver DASSL 
(Differential Algebraic System SoLver ) proposed by Petzold [36] has been found to 
be very robust and easy to use. This method is very suitable for semi-implicit DAE 
(Brenan [37]) such as formed for staged separation process. DASSL is based on back- 
ward difference formula (BDF) which changes its order and step size of integration 
as per the need. The details are described in Brenan [37]. 


(t 


Chapter 2 


Non-Equilibrium Model 
Incorporating Drop Population 
Balance 


2.1 Introduction 


Rigorous simulation of multistage separation processes such as extraction, distilla- 
tion and absorption is, more often than not , based upon the equilibrium approach 
Henley [38]. These models assume that streams leaving stages are in equilibrium with 
each other. Material and energy Balance equations along with constraints such as 
equilibrium relations and normalization equations, are then solved for such models. 
These models are easy to use but predict poor results in practice. 

In actual operation, stages rarely, if ever, are operated at equilibrium despite 
attempts to approach this condition by proper design and choice of operating condi- 
tions. The usual way of dealing with departure from equilibrium is by incorporating a 
stage efficiency into the equilibrium relations. However despite the presence of several 
definitions of efficiencies, as yet no consensus has been arrived at to prove which is 
best. 

In the present work, use of empirical correction factors such as efficiencies or 
HETP is completely avoided; and a non-equilibrium model has been developed for 
the simulation of a counter-current multicomponent extraction process. Since dy- 
namics and mass transfer in a liquid-liquid dispersion are correctly determined by the 
behavior of the dispersed phase, the change in the characteristics (holdup, mean drop 
size d %2 . . . ) of the drop population along the column have to be considered in order 
to describe conveniently the hydrodynamics of column. 

Continuous contact column such as extraction column can be simulated by two 
ways : 
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1. The column hydrodynamics can be represented by plug flow model in which 
dispersed phase (droplets) is considered as plugs of phase moving along the 
height. (Polydispersed model, Jiricny [19]). 

2. Model relying on a drop population balance that leads to the drop size distri- 
bution at every level of the column. Such works have been proposed by Khani 
and Cassamatta [8] for dynamic and steady state simulation of liquid-liquid ex- 
traction (llx) column and Tsouris and Tavlarides for describing dispersion [17], 
holdup profile dynamics and control [16]. 

The above two approaches have been coupled by developing a detailed math- 
ematical model for steady state by Zimmermann [7]. In this model, mass- transfer is 
described by adopting rate based approach such as proposed earlier by Spencer [3] and 
Ricker [2]. A rigorous non-equilibrium stage model as developed by Krishnamurthy 
and Taylor [5, 6] (for distillation column), which has been adopted by Zimmermann 
and his co-workers. 

In present work the basic mathematical formulation of Zimmermann 1992 [7] has 
been adopted to simulate extraction column, both for dynamics and steady state both. 
Drop breakage and coalescence term 11^ is modified from the work by Tsouris [16, 17] 
to facilitate mass transfer. In Tsouris work it was used in drop population balance 
equations for the control of extraction column (such as holdup etc.). Dynamic simu- 
lation has been worked out by including accumulation term in the balance equations 
of above approach. 

Mathematical formulation for steady state simulation is divided into following 
four sections : 

• FLUID DYNAMICS. 

• TRANSPORT RELATIONS. 

• INTERPHASE MODEL. 

• NORMALIZATION EQUATION. 

The model is based on following major assumptions : 

1. Extraction column is operating isothermally. 

2. Stages are assumed to be in mechanical equilibrium i.e. pressure of both phases 
in a stage are equal ( p c = p d )- 

3. Trays are numbered in the direction of dispersed phase flow. This is for conve- 
nience in representation of variables. Variables are defined accordingly so that 
exactly same model formulation can be used whether the dispersed phase is 
heavier or lighter with respect to continuous phase. 
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4. Mixture molar densities are assumed uniform for all drops in the same drop 
class j at any stage k. 

5. No mass transfer resistance at interface between the two phases. 

2.2 Fluid Dynamics 


As with any other model of chemical process, the analysis of the non-equilibrium 
stage starts with the construction of the material balance. Here process is considered 
isothermal 1 . The model as proposed by Khani and Cassamatta [8] is discretized on 
the basis of real stages along column height. 


2.2.1 Conservation relationships 

A schematic representation of non-equilibrium stage is shown in Figures 2.1(a,b). 
Both figures are similar except the flow direction of both phase is reversed. On 
each stage two phases are brought into contact from adjacent stages. The whole drop 
spectrum is discretized in different drop class say 1,2,. . . ,ND based on their drop 
volumes (see section 2.8). Drops of same volume (diameters) are considered as one 
drop class moving as a plug in contact with bulk continuous phase, and the interface 
area is the area between bulk and the drops of that drop class, shown by wavy line in 
Figures 2.1. So in each stage the corresponding two phases are separated by interface 
- thin film of corresponding phases on either side of interface. These films offer 
resistance for mass transfer. Continuous phase is also described as plug flow 2 . Each 
drop class is described by basic variable P^which is the volume fraction of drop class 
j. Summation over all drop classes at any stage k, yields the dynamic holdup (referred 
to as holdup) (j> k of the dispersed phase at any stage k. 

ND 

= ( 2 - 1 ) 
3 = 1 

The mass balance at any stage k for any drop class j and component (species) i 
are as follows : 

MBDijk : Mass balance for Dispersed phase 

~ ( Vfk + Vjiz) x ijk + Vfk+iXijk+i + N ijk + Fj k x{ k + EL ijk = 0 (2.2) 

1 The model can be applied for a column in which at least every stage is considered at uniform 
temperature and temperature gradient along the column may be considered without energy balance 
2 which is interacting with ND plugs of drop classes 




Figure 2.1: Schematic representation of a non-equilibrium stage: (a) Light dispersed 
phase, (b) Heavy dispersed phase. 
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MBCik : Mass balance for Continuous phase 

ND 

Uk-iVik-i ~ (Ujf + U k )yik 4- U k+l y ik+ i — Nij k + Ffy{ k = 0 (2-3) 

j = i 

where subscripts i,j,k takes on value 
i=l,2,...,NC; j=l,2,. . . ,ND; k=l,2,. . . ,NS; 

N^k represents the interface mass transfer rate (moles/s) in a stage from contin- 
uous to dispersed phase, N{j k will be negative if transfer takes place from dispersed 
to continuous phase. is a generation term due to drop breakage and coalescence, 
which is described in detail in section 2.8. 


2.2.2 Flow relationships 

Flow streams V+ k yf k flJ k and Ufl departing from k th stage to the adjacent stage are 
described as a combination of convective and diffusive contribution of flows. Main flow 
direction will contain both convective and dispersion flow; however only dispersion 
flow will exist against the main flow, causing backmixing in extraction column. Axial 
dispersion coefficient D c fc and D d fc are used to represent axial dispersion in continuous 
and dispersed phase respectively at any stage k. 

Flow VPjk : Forward flow for Dispersed phase 

^ = (“* + S ) p ^ s=0 (2 ' 4) 

Flow VMjk : Backward flow for Dispersed phase 

^ = ( ife ) p ^ s=0 (2 - 5) 

Flow UPk : Backward flow for Continuous phase 

u t = (1 - S = 0 (2,6) 

Flow UMk ■ Forward flow for Continuous phase ; 

u t = (1 - M4 S = 0 (2.7) _ 
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Mixture molar densities c k and Cj k are assumed uniform in k th stage for continuous 
bulk phase and all drops of same drop class for dispersed phase respectively. Convec- 
tive contributions are described by the velocity of two phases u k and Uj k respectively 
and diffusive contribution are defined by axial dispersion coefficient. Mixture molar 
densities are given by following relationships. 


d ■- = £l 


(2.8) 

(2.9) 


Velocity of dispersed phase 
ship. 


is expressed by the following kinematic relation- 


4 = -4 + 4 ( 2 - 10 ) 

4 = (i -« m 4 ( 2 . 11 ) 

where Uj k is a relative velocity and u T * k is a characteristic velocity usually available 
as an empirical correlation ( - Laddha [1]). Exponent m ( usually 0-1 ) is a function 

of the hydrodynamics condition, particularly of the Re relative to the drops. In 
present work this is assumed as model parameter and kept unity. However this must 
be evaluated independently by carrying out optimization of model parameter against 
real data for any specific column. 


2.2.3 Transport relationships 

A more rigorous approach is adopted here to represent transport relationships (see 
section 2.9 for details). 

MBICijk : Interphase flow balance, based on continuous phase 

NC NC - 1 

H ij k = Vik ^ 'j Nijk "h aj k ^ ) K^ m j k (jj mk Vmjk) (2.12) 

i — 1 771= 1 

MBIDijk : Interphase flow balance, based on dispersed phase 

NC NC - 1 

Hijk = ~^ijk ^ ' Hijk "b ajk 'y ) Ki m j k (Xmjk %mjk) 

i=l m= 1 


(2.13) 
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Ki m j k and Kf m j k are mass transfer coefficients due to interaction effect in multi- 
component mixture. Evaluation of these coefficients is described in section 2.9.1. 

It should be noted that only (NC-1) independent diffusion fluxes (Taylor [21]) 
exist, since the diffusion flux J* as well as composition gradient must sum to zero 
over the NC species. Interfacial area ajk is directly given by the drop size distribution 
at each stage of the column : 

a jk = 6.0 P jk (2.14) 


2.3 Interface model 


In present model equilibrium between two phases is assumed at the interface with no 
resistance to mass transfer. Equilibrium relation relates the interfacial concentration 
due to both phases by following relation: 


Vijk = Kijkxljk (2-15) 

values of Kij k ■ distribution ratio are calculated by means of ratio between the ac- 
tivity coefficients of both phases. In general Kij k depends upon pressure , temperature 
and composition of both phases at interface. Various standard and reliable methods 
based on group contribution such as UNIFAC, UNIQUAC, NRTL methods can be 
used for evaluating activity coefficients of phases. In present work the simulation is 
carried out with the UNIFAC model (Fredenslund [23]). Various references are also 
mentioned in Bibliography about UNIFAC method and its applicability. 


2.4 Normalization relations 


In order to complete the set of equations corresponding to independent variables, 
normalization equations are written for each phase: 


NDI jk : 


NCI jk : 


NC 



i = 1 


= 1.0 


(2.16) 


NC 

E4* = 1 -° 


i=l 


(2.17) 
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ND jk : 

NC 



^ ^ %ijk = 1*0 

(2.18) 


i=l 


NC k : 

NC 




(2.19) 


2=1 


where (j=l,2,. . 

. ,ND) and (k=l,2,. . . ,NS) 



2.5 Overall simulation model 


Section 2.2 to 2.4 described the model for any general multistage counter-current 
extraction column which is summarized in Table 2.1. The model is discretized in 
such a way that it can be applied for any type of extraction column (agitated and 
non- agitated) such as sieve tray, pulsed sieve tray, RDC etc., without any major 
modification. However correlations used in model may be different for diffrent types 
of columns. In present work the model is analyzed for unagitated sieve tray column 
(counter-current). In such a column there exist two types of holdup - dynamic holdup 
and static holdup. Static holdup has no major relevance in mass transfer except that 
its surface in contact with bulk continuous phase provides mass transfer area between 
the two phases. Further, static holdup causes retardation in the flow of dispersed 
phase, which can be accounted for by velocity correlation of the dispersed phase such 
as proposed by Khani and Cassamatta [8] with the introduction of K v factor in u r * 
correlation. In case of unagitated sieve tray column only dynamic holdup will be 
considered and height caused by the static holdup (beneath or above each sieve tray, 
as the case may be, see Figure 2.2) is deducted from the total compartment height 
of the actual column. The difference is then considered as the active height for mass 
transfer and will be accounted in the volume of stage (14) used for model. Further 
correction in the height may be introduced by the fact that interface area between the 
static holdup and bulk continuous phase also provides area for mass transfer. For this 
area extra height may be added up and used in the volume of stage 14 calculation. 



Feed with solute 


Extract 



Figure 2.2: Schematic representation of extraction column (Light dispersed phase) 
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Table 2.1: Simulation model summary 

^'jk—l^'ijk—1 tyjk T Ijffc )%ijk H” Vjk+\%ijk+l T ^ijk~^~ -^jk^ik 4" Hiifc = 0 

qtiite-1 - (Ut + Ui)Vik + SpViM - Sfafow + FtuL = 0 

v 3 - «+ afe)^*s = o 

*5- (fe)^4s = o 

ttf- (J&Hi-^tes = o 

ur - K+fe)(M0<=|s = o 

Nijk - ~ Oj*Em= 1 1 ^imjk i.Vmk “ ylnjk) = 0 

“ x ijk^2i=l^ijk ~ a jkYLnv=\ K imjk( X mjk ~ x mjk) = 0 
Vijk ~ ■^ijk^'ijk ® 

Effiyfo - i-0 = o 

E^?ary fc - 1.0 = 0 

S&li i-0 = o 

where (j=l,2,. .. ,ND) and (k=l,2,. .. ,NS) 

2.6 Degree of freedom analysis 

Table 2.1 lists the set of equations for single stage, similar sets of equations can also 
be written for other stages, which constitutes the complete problem to be solved. 


For any stage k, following independent variables exist 


Independent variables 

Total No. of variables 

Flow rates of each drop class (V^,V^) 
Flow rates of continuous phase (U k ,U k ) 
Bulk mole fractions (xijk,Vik) 

Interface mole fractions (xfj k ,y(j k ) 

Mass transfer rate (iVy*) 

Drop class volume fractions ( Pj k ) 

Velocity of continuous phase (u%) 

2ND 

2 

(NDxNC)+NC 

2(NDxNC) 

(NDxNC) 

ND 

1 

Total no. of variables/stage 

4(NDxNC)+3ND+NC+3 
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Similarly for any stage following are the number of equations : 


Independent equations 

Total No. of equations 

Material balance equation MBDij k 

(NDxNC) 

Material balance equation MBCik 

NC 

Flow Equations VPj k , VMj k , UPt, UM k 

2(ND+1) 

Transport relation MBIDijk,MBICij k 

2ND(NC-1) 

Equilibrium relations EQijk 

(NDxNC) 

Normalization relations NDIj k , NCIjk 

2ND 

Normalization relations NDj k ,NC k 

ND+1 

Total no. of equations/stage 

4(NDxNC)+3ND+NC+3 


The above problem becomes a set of non-linear algebraic equation for steady state 
ans a set of Differential- Algebraic equations for dynamic state. , which are to be solved 
by any non-linear equation solver ,such as multivariable Newton-Raphson method. 
Solution strategy is described in chapter 3 each dynamic and steady state simulation. 
Providing reasonable initial guess, this method is very robust and converge well even 
for very large equation system. 


2.7 Dynamic model 

Dynamic model is formulated by adding accumulation term in the material balance 
equation sets of steady state model. All other relations in section 2.5 are constraints, 
so they remain same for the dynamic model formulation. Extension of material 
balance sets are given below. 

• DMBDijk 


%-{V k P jk c‘ jk x ijt \ = - O' jt + V£)z + Vf M x ijt+ i + N ijk + Ff t 4 + TUjk 

' DY d 

( 2 . 20 ) 

where DY d (left hand side) term denotes molar accumulation of component i of drop 
class j at stage fc.This term is further simplified as : 


-14 




p jk c%j t x ijk - {M ’ av y k ^[M Ut x iljk j +Cjk x ijkdt P jk 


+ SS d = 0 (2.21) 


• DMBCjk 
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$k) c k Vik] — Uk'-iUik-i {U k + U k )y ik + U k+l yn c + 1 Nij k + F k y{ k 

' v ' i =1 

Dy - s 

( 2 . 22 ) 

where DY C (left hand side) term denotes molar accumulation of component i at stage 
k in continuous phase. This term is further simplified as : 

-U [d - I + I jM +ss ‘ = ° 

(2.23) 


All other equation sets described in section 2.5 will remain same for dynamic 
modeling also. These relations act like constraints at any instant of time. Degree of 
freedom analysis discussed in section 2.6 is still valid since number of equations and 
number of independent variables remained unaltered. 


2.8 Drop Breakage and Coalescence 

Smooth flow of dispersed phase in form of droplets in extraction column can be de- 
scribed by the simple transport relations. Abrupt changes in drop distribution due 
to drop breakage and coalescence cannot be handled by the transport relations alone 
satisfactorily, instead population balance approach is introduced for such situations. 
Especially in case of extraction where large transport area is created by agitation, 
accompanied by high drop breakage and coalescence, need of such model is immense. 
Because of breakage and coalescence phenomena dispersed phase droplets continually 
forgo their individual identities. The material domain for which the transport equa- 
tions are applicable is thus subject to discontinuous changes in an apparently random 
manner. 

Ramakrishna [9] gives an excellent review of the several works that have been 
published on drop breakage and coalescence during past 2-3 decades. Among the first 
published works in the field, Coulaloglou [10] proposed a phenomenological model to 
describe drop breakage and coalescence in a turbulently agitated liquid-liquid disper- 
sion. It was a theoretical model based on fundamental concepts. The same model has 
been worked out by several other workers Sovava [11, 12], Rubio [15], Hsia and Tavlar- 
ides [13, 14] and Tsouris [16, 17]. These works utilized the drop population balance 
ignoring mass transfer aspects. Such a model has been well applied by Tsouris [17] 
for drop-breakage analysis, holdup profile and holdup control of the column. The 
work presented here incorporates the concept of these models (based on number of 
drops balance) into the mass transfer model in this work. Drop breakage rate II-* 
and coalescence rate Ilgf terms are modified as described in the next paragraph. 

Earlier drop breakage and coalescence rates were written for population balance in 
which balance over number of drops was used. In present model number of drops and 
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their concentration (of each drop class) are easily calculable as simulation variables. 
So these rate terms in form of number of drops per time as described in Tsouris [16, 17] 
are simply multiplied by the volume, molar density and concentration of respective 
drop (drop class) to get corresponding mass transfer rate (moles/s) due to drop break- 
age and coalescence. These transfer terms in all are called as generation term Ily* 
due to drop breakage ngf and coalescence ngf and thus included in material balance 
equation set MBDijk see section 2.2. However for non-agitated columns generation 
term may not be very much significant at any stage of extraction column. 

Following sections describe the dichotomy of Uij k term and its derivation from 
population balance model. Various correlations were also used in the original works, 
some of them still need detailed research, in fact drop breakage and coalescence is 
still a topic of intense research at present in the field of liquid-liquid dispersion. 

2.8.1 Generation due to drop-breakage and coalescence: 

In dispersion drop spectrum is discretized on the basis of volume of drops (say 0 to 
Vmax ) and each discretized part is assigned a number as 1,2,. . . ,ND. Diameters are 
than calculated for each drop class. The generation term n^* (negative for loss) in 
material balance equation MBDijk is divided in two parts as below. DB and DC 
refer to drop breakage and drop coalescence respectively. 

n«i, = ngf + ngf (2.24) 

Each of two terms on the right hand side of above equation is further divided 
into two parts as actual loss and actual gain of material from any drop class at any 
stage. Subscripts i,j,k are used for component, drop class and stage respectively. 
For all terms base formulation are referred from Tsouris [16], Coulaloglou [10] and 
Sovova [11]. They are modified for mass transfer here. 

Generation due to drop breakage : ngf 

In present formulation sometimes only drop class is referred (by j), but in actual 
calculation either drop volume or drop diameter is used whichever is convenient. The 
term ngf is divided as : 

ngf = -ngf + ngf (2.25) 

• ngf" ( Material lost from jth drop class due to breakage of drops of jth class ) 
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This material will be distributed in lower 3 drop class. For lowest drop class this 
term is not applicable. 

= [9kti)nkti,t)]vjCj k x ijk (2.26) 

Correlations for gk(j) (breakage frequency) can be referred from Appendix B. 
n k(j,t) is a number of drops in drop class j, at any time t. Vj is a volume of 
single drop of drop class j . 

• Ilg£ + ( Material gained by jth drop class due to breakage of higher i drop class 


This material is gained by the drop classes (except highest drop class) due to 
breakage of other drops. 

_ . rvmax 

ngf = f [Pivj^'^^giv^niv'.t^Vjcj^Xijkdv' (2.27) 

J Vj 


Breakage type u is assumed as binary breakage for present study. A daughter 
probability density function P(j,j') is used to obtain probability of breakage for 
drops of volume v into Vj , which must satisfy some conditions, Ramkrihshna [9, 
page 58] such as 


P{di,dj ) = 0.0, 



ifi = j{ordj < di ) 


(2.28) 

(2.29) 


where j= father drop class i= daughter drop class 


Appendix B can be seen for 0(vj,v') correlations. 

Generation due to drop coalescence : Ilgg 
Similarly this term is divided into two parts. 

ng? = -ngf + ngf (2.30) 

• ngg - ( Material lost from jth drop class due to its coalescence with other drops) 


3 means drop class having less diameter of jth drop class 
4 means, drop class having large diameter than jth drop class 
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This material will go to other possible higher drop classes. 

j (“Vmax~~Vj 

ngj? = n(v j ,t)v j c d jk x ij k X (v^v^hiy^v^n^dv' (2.31) 

Vmin 

Since after coalescence drop volume is not expected greater than v max hence 
upper limit of integration is set as (v max — Vj). 

• ngf + ( Material gained by jth drop class due to coalescence of other small drops ) 


This material is gained by the drop class (except lowest drop class), due to 
coalescence of other smaller drops. 


ngf f = f 3 [A(u y - v', v')h(vj - v\ v')n{vj - v', t)n{v ', t)] Vj(c%)' {x ijk )'di J 

J Vmin 

(2.32) 

where drop molar density of the new drop formed by coalescence 


/ <t y fjfck] v> + ~ ^ ; ] Z 

v' + (Vj — V 1 ) 


and its molar fraction 


i X ijk) 


u jk 


V 


' Xijk[v'] + C d jk [Vj - V 1 ] (vj - v') x ljk 


Vd — V 7 


CjkW] V + C d jk [v j - V'} {Vj - v') 


(2.33) 


(2.34) 


NOTE : Drop volume enclosed in brackets [ ] in abo ve equ ations represents the 
corresponding drop class for that variable. For example c d k [v'} represents the molar 
density of drop class corresponding to drop volume v' . 


h(vi,Vj) is the collision frequency of two drops and A(uj. Vj) is the coalescence 
efficiency between two drops. The efficiency \(vi,Vj) term is not well understood in 
liquid-liquid dispersions and various correlations are proposed by different authors. 
Some favor high efficiency for larger drops and some for smaller drops. h(yi , Vj ) term, 
originally derived from kinetic thery of gases (Kennard [39, pp 103-104], has been 
further modified by other authors. These correlations are listed in Appendix B. 

The integration terms in above formulation are evaluated by using trapezoidal 
rule, since discrete points are already known in form of distinct drop class (drop 
diameter or volume) from 0 to v maXl hence no separate discretization is necessary for 
trapezoidal rule. Similar formulation based on trapezoidal rule for drop population 
model can also be seen from Sovova [11], which also described explained discretization 
errors. 

In general drop generation term Ily*, is function of : 


(2.35) 
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where bold subscript j represents the general functionality valid over j=l,2,. . . ,ND 
, whereas i,k are same in both right and left hand side of above functionality relation. 

2.8.2 Discretization and drop size distribution at inlet 

For simulation either a single drop class is assumed with average drop diameter (sauter 
mean diameter) in dispersion calculated by suitable correlation (Laddha [1], Baird [4]) 
or a discretized drop spectrum is taken. For dynamic simulation a single drop class 
approach is advisable computationally. However for a broader analysis of extraction, 
one needs to include more than one drop class in simulation. More number of drop 
classes means less discretization error but at the cost of higher computational time. 
Optimal number of drop classes are used for simulation based on these two factors. 
Usually about 8 drop classes are enough to avoid any significant discretization error. 

As described earlier discretization is made on the basis of drop volume usually 
between v m i n (= 0) to v max ; where v max is calculated by a suitable correlation such 
as Baird [4, 1983, pp 129, eqn 17]. Let’s say v min and To are the minimum and mean 
volume of drops respectively in a given dispersion, then a discrete point for n th drop 
class is given by: 

v(n) = v min + n* Av (2.36) 

The feed distribution can be assumed as a Gaussian probability density distribu- 
tion function which suggests that the probability of drops falling with in range v to 
v + Av is given by. 

Aw= v^ exp H(^) 2 ] (2 - 37) 

where — oo < v < oo 

here a is the standard deviation. 

So the number fraction of drops that may lie from (j-1) to j th drop class but 
accounted in j th drop class is given by 

Pj = f A(vj)dvj (2.38) 

J Vj-l 

which must satisfy for real problem, 

rvmax 

/ A(vj)dvj = 1.0 (2.39) 

Jo 

if above is not satisfied then normalize to give. 

_ j:u 

Pi JT“ 


(2.40) 
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Figure 2.3: Normal distribution for inlet drop spectrum (Volume probabilty density 
vs volume of drops) 

N 0 the total number drops at inlet is calulated on the basis of mean drop volume, 

JV„ = (2.41) 

^0 


here f o the holdup is assumed for initial calculations. Volume fraction of drop 
class j for inlet is calculated as : 


PjO — 


{NoPj)vj 

Vi 


(2.42) 


substituting No, we get 


PjO = 


<?0PjVj 

vo 


(2.43) 


In order to avoid any anomaly due to discritization this can be normalized again 
for getting Zf=o Pjo = io 

Further inlet distribution of dispersed phase flow rate is decided as 




V 


. Pr 


jo 


<f>0 


(2.44) 


for all j=l,2,. . . ,ND 

This way discretization and inlet distribution of flows are fixed. Usually inlet 
distribution depends upon the way inlet is introduced to column, i.e. by using nozzle, 
sparger etc. If somehow, actual distribution is known that can be used instead of 
theoretical distribution suggested above. For other stages 1 to NS , the variable Pjk 
is calculated and simulation results are used to find out drop distribution at any stage. 
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2.9 Estimation of Mass Transfer Coefficient Ma- 
trix 


Recent developments and success of rigorous simulation model, using coupling be- 
tween various species present in separation processes such as distillation and absorp- 
tion Krishnamurthy [5, 6], encourage the use of such rigour in modeling of extraction 
also. 

Mass transfer in multicomponent mixture is more complicated than in binary 
system because of possible coupling between the individual concentration gradients. 
Phenomena such as reverse diffusion 5 or osmotic diffusion 6 are possible in multicom- 
ponent system (Taylor [21]). Models of mass transfer that are able to account for 
interaction effect are described here. 

Mass transfer is a rate process driven by gradients in concentration. The rate 
of mass transfer in the continuous and the dispersed phase depends upon the con- 
centration gradient across the film. Here bulk concentrations are assumed uniform 
and equilibrium is assumed at interface. Many other factors such as temperature, 
concentration, physical properties of the fluids (density, density difference, viscosity, 
interfacial tension etc.), system geometry and hydrodynamics , influence the rates 
that are achieved. Current practice is to lump the effects of most of these variables 
into mass transfer coefficient. Mass transfer in a multicomponent mixture is descibed 
[Ni = yi J2i> Ni' + J*] as in Bird [40] and Taylor [21]. On right hand side of this expres- 
sion, first term represents the convective mass transfer and second term represents 
the diffusive mass transfer. 


2.9.1 Estimation approach 

A more rigorous approach is used to write transport relationships MBIDijk and 
MBICijk ■ Resistance to mass transfer (offered by both phases) can be accounted for 
by using separate rate equation incorporating interactive mass transfer coefficients. 
Interfacial area between bulk continuous phase and any drop class is easily calculated 
by equation 2.14 as a^. 

Methods that take interaction effects into account, but are implicit in fluxes or 
transfer rate 7 , are proposed by Krishnamurthy [5, 6] for distillation and absorption. 
Same philosophy is used here, but variables are tailored for extraction. In general, 
transfer rate for continuous phase is written as (similar development is true for dis- 
persed phase): 

5 diffusion of species against its own concentration gradient 

6 diffusion of species even though no concentration gradient for that species exist 

7 Here transfer rate as Nijkis used directly, since in material balance actual transfer rate is used 
instead of flux 
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nc ; 

Hijk — TJik ^ ' N i'jk “I" &jk J ijk (2.45) 

i'=l 


Flux J* jk can be written by incorporating interaction effect into account, by us- 
ing generalized Fick’s law, Taylor [21], based on mass transfer coefficient instead of 
diffusivity. Equation 2.45 can be rewritten as : 

NC NC - 1 

Nijk = Uik 22 Ni'jk + a jk 22 ^imjk^Vmk ~ Umjk ) (2.46) 

i '= 1 m= 1 

In matrix form above equation can be written for each drop class “j” and in stage 
k th as: 


(N)jk = (y)jkNt + a jk [K c ] jk (y - y ! )jk (2.47) 

Where Nt = N Vjk 

() represents vector 

[] represents matrix in equation 2.47 

Equation 2.47 consists of (NC-1) equations for any drop class j at stage k , since 
Eila J*jk = 0-0 • hence only (NC-1) independent equation will exist. 

In general mass transfer coefficient k is dependent on 
k = {(bulk cone., Interfacial cone., flux or transfer rate ) 

. Mass transfer coefficient matrix of size (NC-1) by (NC-1) are then evaluated as : 

[K c ] jk .a jk = [B]-£a jk [$] jk {ex p[$] ifc - [J]}” 1 (2.48) 

This expression is then evaluated by following one of the procedures. 

• Krishna and Standart method This method is used in present study. 

= [M]j k with m ijk = y ik (2.49) 

ajk 

=: with TTiijk = iV ijk (2.50) 

• In linearized theory due to Toor(1964) 

\B ] i 1 

i- = [JVf] jk with Triijk = yik = ~( Vik t Vink) 

. ■ CLjk ■ ■ ■ " . 


(2.51) 
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Where general matrix [M]j k components are evaluated as given below : 


Maj k 


m ijk 


M-im'jk — 


-mijk 


k-iNCjk a jk 

1 


NC 

+ £ 

m / =l 

mJ^-i 

l 


YYl m f jk 
_. x k im/jk 


v kira'jA; k k^ NCj k ®jk , 


z = l,2,...,(AC-l) 

i = l,2,...,{NC-l) 


(2.52) 

(2.53) 


Above method is described in full detail by Taylor [21, Chap. 7,8] and is based on 
film theory. General solution of Krishna and Standartis is applicable to mixture with 
any number of species. Some limiting cases are also described there. Zimmermann [7] 
used linearized theory of Toor(1964), whereas Krishna and Standart is used here. This 
method is general and gives exact solution. 


2.9.2 Evaluation of binary pair k im jk from correlation 


The binary local mass transfer coefficient k im jk for any drop class j and stage fc, 
can be calculated by various empirical correlations (see Appendix B) available in 
literature [4], Rousseau [30],Treybal [41]. Usually these correlations are available in 
form of Sherwood number (Sh) relationship, such as : 

Sh = ci Re C2 Sc c 3 (2.54) 

Where Ci,C 2 ,C 3 are empirical constants. 

Sh = (2.55) 

Pm 

where = binary mass trans fer coef ficient 
p m — molar density = Cj k or c k 

Vi m — Diffusion coef ficient [Maxwell) of binary pair i — m 

In multicomponent mixture diffusion coefficient Vi m is given by (see equation 
4.2.18;Taylor [2f.:pp 91, is rewritten here) 


V ir , 


(Pin)' 


W«*Y 


(2.56) 


Where infinite dilution diffusivity V° is evaluated by Wilke - Chang (1955) corre- 
lation see Taylor '21, pp 73] or Reid & Sherwood [42, pp 549]. Only mutual diffusivity 
V° m (i ^ m) needs to be evaluated here. However for self diffusion coefficient sepa- 
rate correlations are used, though these are not required in section 2.9.1 calculations. 
Binary mass transfer coefficients calculated here are used in section 2.9.1. 



Chapter 3 


Steady State and Dynamic 
Simulation 


3.1 General 

Section 2.5 desribes the overall simulation model to be solved for independent sim- 
ulation variables listed in section 2.6. The simulation model described is a set of 
non-linear equations which has to be solved by using multivariable Newton-Raphson 
method. A computer code is written in “C” langauage for this purpose, both to 
solve steady state and dynamic simulation of multicomponent unagitated sieve tray 
column. However with the replacement of correlations, the same model as well as com- 
puter program can be used for other types of extraction columns (agitated columns) 
also. Newton-Raphson method for solving set of non-linear equation is found to be 
very stable. Extension of steady state simulation is used for dynamic simulation by 
using DASSL (Differential Algebraic System SoLver). In dynamic simulation some 
elements (related to balance equations) of Jacobian is evaluated by using DASSL and 
all Jacobian elements used for steady state are used without alteration. The system 
is again solved using Newton-Raphson(N-R) method at any instant of time. 


3.2 Steady state simulation 

Let (f S ' S {X)) denotes the a vector of non-linear equations to be solved ( as simulation 
model listed in Table 2.1) for collection of independent simulation variables X = 
(J s f 1 ,A' 2r ...,A' n ), with solution (f s _ s (X)) = (0) . Then N-R method suggests the 
solution improvement X p+l by the following relation. 


- X-) = -(U*(X r )) 


(3.1) 
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Where [J p ] s . s is the Jacobian matrix with elements : 


r Bfj (X) 

dXj 


at any iteration “p" 


(3.2) 


where Jij G [J p ] s .s and fi(X) G (f s . s (X)) . This method is very robust and versatile 
except that it requires the calculation of the Jacobian which involves derivatives of 
equations to be solved with respect to independent variables. These independent 
variables are listed in section 2.6 and here denoted by single vector X. 


3.2.1 Jacobian evaluation for steady state simulation 

The analytical derivatives are evaluated by partial differentiating the model equations 
(as in Table 2.1) with respect to independent variables as denoted by X. Derivatives 
of some complicated terms (such as such as ^^(distribution ratio), mass transfer 
coefficient Kf m j k ,Kf m: j k , drop breakage and coalescence generation term ITy*) with 
respect to independent variables (X), which are difficult to derive analytically are 
evaluated numerically. Since calculation of numerical derivative is a time consuming 
process, good bookkeeping is necessary during simulation to avoid duplicate calcula- 
tions and to save simulation time. In order to improve simulation time further, some 
other hybrid techniques are proposed ( Krishnamurthy [5, 6]) in literature. Sometimes 
however these methods may cause an instability problem due to less accurate update 
of Jacobian elements. Since in this work accuracy is stressed more, derivatives are 
evaluated analytically to the extent one can evaluate comfortably and the remaining 
derivatives are evaluated numerically. For storage and calculation variables double 
precision arithmetic is used on Unix based machines (DEC a 2000, HP9000-850 m/c). 

As for any other staged separation process, a Block-Tridiagonal Jacobian matrix 
(see Figure 3.1) is formed. Where block “A” represents the interaction from previous 
stage, block “C” represents the interaction due to to next stage and block “B” is the 
main block for any stage in Jacobian [J s „ ,]. Apart from blocks A,B and C remaining 
Jacobian matrix is full of zeros. Vector D contains the residuals of the equations to 
be solved by Newton-Raphson method. 

The extent of sparsity increases with number of stages. Hence a linearized sim- 
ulation model equation 3.1 is solved by some suitable linear equation solver such 
as Block-Thomas algorithm [43]. Block-Thomas algorithm requires only storage of 
blocks A.B,C and function vector D for every stage. The size of all A,B and C blocks 
is same as the square matrix (m x m ) , where m=number of independent simulation 
variables per stage (see section 2.5). Residual vector D will be having size of (m x 1). 

In order to make evaluation and coding of computer program easier each block 
(A,B or C) is further divided into subblocks or minors (M^). Selection of rows and 
columns of Jacobian for each minor is shown in Figure 3.2. The rows of minors 
are marked (by the type of equation sets, as described in section 2.5 and column of 
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Figure 3.1: Block Tridiagonal Jacobian matrix (superscript represents stage number) 
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minor matrix is marked by class of independent variables described by section 2.6. 
A general block containing 30 minors is then represented as in Figure 3.2. These 
minors contains partial derivative of equations (listed as row heads) with respect to 
independent variables (listed as column heads). 

As per “C” language convention minors subscripts are starts from (0,0) intead of 
(1,1). Similarly vector D is divided into 5 (5 x 1) parts. These minors and vector D 
are then calculated separately for one stage and finally synthesized as A,B,C blocks 
and D vector respectively for that stage. 

Steady state simulation strategy is then summarized as below : 

NOTE : Stages are numbered along the direction of dispersed phase flow. 


# 1 : Read data (Computer file : simulfun02.c (Table C.l) 

— » Column specific (design parameters) as in Table C.2 

— ► System specific (Physical properties, molecular weight etc. as in Table C.3) 
— *• Operating condition (P,T,flow etc. as in Table C.4) — + Read UNIFAC pa- 
rameters as in Table C. 

# 2 : Condition input data (simulfun02.c Table C.l) as per the simulation variables 

used for column simulation, i.e. feeds are represented in terms of variables 
used in simulator. Carry out discretization of drop spectrum based on the 
drop volumes. This discretization (ND) and number of stages (NS) along with 
number of components (NC) decides the total number of independent equations 
and variables to be solved. 

# 3 : Set iteration number, p «— 0 . 

Suggest initial condition as : 

For main flows ( V + ,U~ ) — » Take same as total feed flow up to that stage for 
respective phase. 

( Xijk,Vik ) — ► same as feed concentration. 

(xlj k ,ylj k ) — ► same as feed concentration. 

(Nijk) — » very low value such as 10~ 06 . 

(Pj k ) — ► Assume some value between 5-15% 

(u%) — > depends upon column and flows , typical value is (3mm/sec). 


# 4 : Increment p +— p+1. 

Solve set of simulation equations (step.c Table C.l) for independent simulation 
variables by Newton-Raphson method, at step "p” and update simulation vari- 
ables as 


* p+1 = x ,_ A TOi(/„(^)) (3.3) 

Where A is calculated by single variable optimization (Golden section search) 
of objective function (i.e. sum of square of residuals of model equations) . 




(General Block A, B orC) (Corresponding 

part of D vectors ) 
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# 5 : If termination condition satisfied then stop simulation else go to step # 4. 

# 6 : Report solution at steady state. 


3.3 Dynamic simulation 

In order to distinguish between steady state and dynamic simulation formulation “Y” 
is used for independent variable vector at any instant of time in place of “X” used for 
steady state. By modified material balance for dynamic modeling DM BDijk , DMBCik 
(as in section 2.7 the new set of equations can be written as : 


fd.s{tn+ 1) y.+i,K +1 )=0 (3.4) 

were Y’ is a derivative vector containing time derivative of independent variables. 
Equation 3.4 represents the Differential-Algebraic Equation (DAE) system of index 
more than 1 (using terminology of Brenan [37]). 

In order to solve equation 3.4 one can use backward difference (BDF) technique 
such as Euler’s method. For Euler’s method time derivatives Y’ vector elements ( y ! 
) are represented as : 

K +1 = (^) (3.5) 

“h” is current time step and y is the element of vector (Y), thus DAE becomes at 
any instant t n+l 

fd.s (t n + 1, y n +u = 0 (3-6) 

This equation could be solved for Y n+ i at any instant of time t n+1 by using previous 
known value Y n . with starting Yo as initial guess. Other techniques such as Gear’s 
algorithm, Runge-kutta method could also be applied to solve DAE problem. 

Since, steady state model has been solved by N-R method, it was found conve- 
nient to use any BDF technique which could be solved by N-R method with little 
modification. Euler’s method is very inefficient to solve such a problem. A new ro- 
bust BDF technique known as DASSL (Differential Algebraic System SoLver) is used 
as proposed by Petzold [36, 1982]. This method is quite robust for DAE system for 
index 1 or lower and works efficiently for semi-implicit DAE such as present dynamic 
model formulated for LLX column. For using DASSL formulation initial condition Yo 
and Y 0 ' are required. It is assumed that DASSL is applied from one steady state to 
other steady state. Hence, when first DASSL is applied Y 0 ' = (0) is assumed; if not, 
then either these values have to be supplied or calculated (as in Brenan [37]), 

By using DASSL formulation (as given in Appendix A 

7 '( 0 ) a s / ( 0 ) V 

Vn+1 Un + 1 l {.Vn+l Vn— l) 

O'n+X 


(3.7) 



3.3 Dynamic simulation 


31 


This is substituted in equation 3.6 and solved for vector Y n+ \ by using N-R method 
at any instant of time. 

The N-R problem at any instant of time f n+1 is given by: 


[j’V.K+i - vi) = -{/„,(n) 


The elements of new Jacobian [J p ]d. s 

i d J i¥l 


are given by 
at any iteration 


Y 


(3.8) 


(3.9) 


where fi(Y) € fd.s and J l3 € [J p }d.s- The Jacobian for dynamic simulation problem 
is then given as : 


| [J p 


\d.s 


dfd.3 


dY 


dfs.s 

dY 


+ 


df±s 

dY 


+ a 


dhs 

dY' 


(3.10) 


Here independent variable vector (Y) is exactly similar to vector (X), desribed in 
steady state formulation except (Y) is at any instant of time f n _^. Where matrix 
contains derivative of accumulation terms only, which includes derivatives 
of DYd and DY C with respect to elements of vector (Y) only. Derivatives of these 
accumulation terms with respect to time derivative (Y’) are included in 

Derivative of remaining terms are already included by matrix 


df,, 


dY 


dY' 

So, 


matrix. 


[J p \ 


d.s 


df d .s 


dY 


-fa 


dUs 

dY' 


(3.11) 


where 


'df d . s ' 


dfs. s ' 


'dfd.3 

_ dY _ 


. dY _ 

+ 

. dY . 


(3.12) 


For dynamic simulation only elements of and a j have to be evaluated 

extra and will be appended to the formulation of Jacobian for steady state. 


DASSL is a variable order and variable step size algorithm which makes it robust. 
When system tends towards steady state, step size increases automatically and at 
perturbation it reduces the step appropriately. Formulation of DASSL, step size and 
order selection strategy can be seen in Appendix A. 



Chapter 4 


Results and Discussion 


4.1 Application of the model 

In principle the model described in previous chapters can be applied for any type of 
multicomponent, staged countercurrent extraction column. In this thesis the validity 
of such a model is tested for unagitated sieve tray extraction column. Two extraction 
systems have been tested in detail; results of which are presented here. 

System 1 n-Heptane(l),cyl-Hexane(2),Benzene(3),Toluene(4),o-Xylene(5) and Sul- 
folane^). 

This is an industrial extraction system commonly known as aromatics extrac- 
tion. Here sulfolane as a solvent, extracts BTX (Benzene.Toluene, Xylene) from 
naphtha feed, which contains n-heptane, cyclo-hexane along with BTX. Sul- 
folane, as a solvent is kept in dispersed phase and fed from the top (first tray) of 
the 40 staged column see Appendix C; whereas light naphtha feed in continuous 
phase enters from the bottom (40th tray) of the column. Products leaving the 
column are — Raffinate phase (purified naphtha) from the top of the column 
and Extract phase (which contains solvent sulfolane and BTX extracted from 
naphtha) from bottom of the column. In real operation this extract phase is 
further washed with steam, BTX is separated out and purified sulfolane recir- 
culated to the extraction column. Data for a 30 cm diameter and 40 stage 
unagitated sieve tray column in this study have been obtained from Engineers 
India Limited (EIL), Gurgoan. Column design parameters, system properties 
and operating conditions obtained from EIL are listed in Appendix C. These 
data are then fed to the simulator program developed as a part of this thesis 
and results obtained. Naphtha as a continuous phase feed enters with composi- 
tion (0.37.0.37, 0.05,0.16,0.05,0.0) x ; BTX in feed is extracted by pure sulfolane 
in dispersed phase having feed composition (0.0, 0.0, 0.0. 0.0, 0.0, 1.0). Extrac- 


: 6 components, described as System 1 
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tion operation has been considered isothermal at 95°C', thus physical properties 
are evaluated at this temperature. However, UNIFAC parameters available at 
25 ° C in literature are used in the absence of any other reliable data. UNIFAC 
parameters reported by Mukhopadhvaya [28] have been tried out also in this 
work. These parameters have been evaluated at 100°C for some organic com- 
pound, however use of these parameters only favors the extraction of Benzene 
and Toluene and Xylene remained unextracted which is not true for practical 
column consodered in this study. The reason of such failure may be that these 
parameters are evaluated based on binary or tertiary mixtures of limited compo- 
nents whereas data reported in Magnussen [24] (used finally in this work) have 
been evaluated based on binary or tertiary mixtures of thousands of componets. 

System 2 Water(l),Acetone(2),Toluene(3). 

This is a test system recommended for the lab experiments. Toluene as a solvent 
(dispersed phase) enters from bottom of the 4 tray unagitated sieve tray column, 
which extracts, Acetone from the continuous aqueous phase entering from top 
of the column. Data are listed in Appendix D and results are also included 
therein. 


4.2 Steady state results 

Test SYSTEM 1 

For the ease of simulation a single drop class having an average drop diameter is 
considered initially with no breakage and coalescence effects. A total of 1440 (36 
equations per stage) equations are solved by using Newton-Raphson (N-R) method, 
with convergence criteria as sum of squares of residuals of equations being less than the 
given tolerance 10 -10 . N-R method converges normally in 4-5 steps for all simulation 
runs described here. Data for all simulation runs are listed in Appendix C (Table C.2 
to C). Here 5 simulations runs (run 1 to 5) for steady state are presented. Feed 
conditions for each of these simulation run are listed in Table 4.1. 

Simulation run 1 (Figures 4.1 to 4.7) 

Figures 4.1 and 4.2 show the concentration profiles of all 6 components in continu- 
ous and dispersed phase respectively. From theses profiles it is clear that sulfolane 
extracts almost all BTX contents from the naphtha feed. Naphtha product as a raf- 
finate phase leaves with 50.5 % n-heptane and 47.6 % cyl-hexane. Extract phase 
is leaving with 87.4 % sulfolane, remaining being BTX with negligible amounts of 
cyl-hexane and n-heptane. Solutes BTX concentration profiles are further enlarged 
in figures 4.3(a,b). Sulfolane selectivity is highest for Benzene then Toluene and lowest 
for o-Xylene. Hence all benzene is extracted up to 30 th stage (i.e only in 10 stages from 
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Table 4.1: Feed conditions for steady state simulation runs 


Run 

no. 

S 

( moles /s) 

F 

(moles/s) 

S/F 

(w/w) 

Naphta feed concentration 
(mole fraction) 

Remarks | 

1 . 

4.5 

2.0 

2.7 

0.37,0.37,0.05,0.16,0.05,0.0 

Base case 

2. 

4.5 

3.0 

1.8 

0.37,0.37,0.05,0.16,0.05,0.0 

S/F 

decreased 

U 

5.5 

2.0 

3.3 

0.37,0.37,0.05,0.16,0.05,0.0 

S/F 

increased 


6.75 

3.0 

2.7 

0.37,0.37,0.05,0.16,0.05,0.0 

Throughput 

increased 

5. 

4.5 

2.0 

2.7 

0.37,0.37,0.09,0.10,0.07,0.0 

BTX 

cone. 

changed 

10. 

4.5 

2.0 

2.7 

0.37,0.37,0.05,0.16,0.05,0.0 

With two 
values 
of axial 
dispersion 
coeffi- 
cients 

D c = 3.0 
and 

35 cm 2 /s. 


S=Solvent feed rate (Average molecular weight = 120) 
F=Naphtha feed rate (Average molecular weight = 100) 
S/F= Solvent to feed ratio 
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feed inlet), whereas toluene and o-xylene require full column height to be extracted 
fully. 

One of the key variable in the simulation of extraction process is a dynamic 
holdup, which is found in the range of 5.2 % to 5.8 % throughout the column see fig- 
ure 4. 4. Forward (main) flows of both phases are also plotted in figures ??. For present 
case continuous phase axial dispersion coefficient D c = 3.0 cm 2 /s is assumed, whereas 
no axial dispersion between stages is assumed for dispersed phase i.e. D d = 0.0cm 2 /s. 
Effect of this dispersion is clearly seen in figure 4.6, which shows higher forward flow 
in the intermediate stages as compared to inlet and outlet flow rates. Continuous 
phase velocity which is one of the independent variable is also plotted against stage 
number in Figure 4.7, ranging from 2.6 to 3.8 mm/s. 

Profiles shown and described above are in agreement of real aromatic extraction 
process. However exact value may not match at this point of time, since model is not 
fine tuned against real plant data.In the absence of any real plant data, against which 
model could be tested, some logical tests can be performed over simulation package 
to check its validity. Simulation run no. 2,3,4 and 5 are devised for this purpose. 
Simulation runs are performed with the data set corresponding to these runs and each 
of these runs are then compared independently against results obtained in simulation 
run 1 (base case). 


Simulation run 2 (S/F ratio decreased, Figures 4.8,4.12,4.13) 

In this continuous phase inlet flow rate (Uin) is increased from 2.0 moles/s (0.20 
kg/s) to 3.0 moles/s (0.3 kg/s), keeping solvent (dispersed phase) flow rate as con- 
stant at 4.5 moles/s (0.54 kg/s), thus decreasing S/F ratio 2 from 2.7 w/w to 1.8 
w/w. Figures 4.8(a,b), show the concentration profiles of BTX (solutes) in continu- 
ous and dispersed phase respectively. Since S/F ratio is decreased it should cause less 
extraction. This is clear from figure 4.8(a), which shows naphtha product leaves as 
a raffinate with higher aromatic contents (B=0.0%, T=3.1%, X=2.7%) as compared 
to in figure 4.3(a). On the other side extract phase also leaves with higher content 
of extracted aromatics (since feed contents of BTX has increased) as compared to 
figure 4.3(b). Dynamic holdup ranges from 5.3% to 6.0%, (as seen in figure 4.12) not 
much change is expected, because there is no change in the flow rate of dispersed 
phase. However continuous phase velocity registers a change, which peaks up to 6.0 
mm/s as seen in figure 4.13 compared to 3.8 mm/s in Run 1. 


Simulation run 3 (S /F ratio increased) (Figures 4.9,4.12,4.13) 


In this set, dispersed phase feed flow rate is increased from 4.5 moles/s (0.54 kg/s) to 
5.5 moles/s (0.66 kg/s), keeping continuous phase flow rate constant correspond to 


2 Solvent to feed ratio 



MOLE FRACTION MOLE FRACTION 
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Figure 4.1: Concentration profile for dispersed (Sulfolane) phase (Run 1) 



Figure 4.2: Concentration profile for continuous (Naphtha) phase (Run 1) 
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Figure 4.5: Forward flow (continuous phase) (Run 1) 





MOLAR F 






MOLE FRACTION MOLE FRACTION 



Figure 4.8: Solutes (BTX) concentration profiles (a) Continuous naphtha phase (b) 
Dispersed sulfolane phase. (Run 2) S/F decreased 
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simulation run 1 as 2.0 moles/s (0.20 kg/s), thus increasing S/F ration from 2.7 w/w 
to 3.3 w/w. Again Figures 4.9 (a, b) show concentration profiles of solutes (BTX) in 
continuous and dispersed phase respectively. As per expectation increased S/F ratio 
causes early extraction of solutes thus reverting concentration profiles for toluene 
and xylene in continuous phase. Solvent phase even after full extraction carries lower 
percentage (although high amount) of BTX, as compared to in figure 4.3(b). Dynamic 
holdup ranges from 6.4% to 7.1%; it is up by 2.0% approximately, mainly due to 
increase in solvent flow rate as seen in figure 4.12. Continuous phase velocity doesn’t 
register any significant change as seen in figure 4.13. 

Simulation run 4 (S/F constant, throughput increased) (Figures 4.10,4.12,4.13) 


For this run, both continuous and dispersed phase feed flow rates are increased as : 

Continuous phase — from 2.0 moles/s (0.20 kg/s) to 3.0 moles/s (0.30 kg/s). 
Dispersed phase — from 4.5 moles/s (0.54 kg/s) to 6.75 moles/s (0.81 kg/s). 

So flow rates are increased but S/F ratio is kept same at 2.7 w/w. Since column 
throughput is increased, solute concentration profiles do not change as is clear from 
comparison of Figures 4.10(a,b) to Figures 4.3(a,b). However both dynamic holdup 
and continuous phase velocity register an upward change in their values. Dynamic 
holdup (Figure 4.12) peaks to 9.2% and velocity to 6 mm/s (Figure 4.12). 

Simulation run 5 (BTX concentration changed) (Figures 4.11,4.12,4.13) 

In this set, concentration of naphtha feed is changed from (0.37,0.37,0.05,0.16,0.05,0.0) 
to (0.37,0.37,0.09.0.10,0.07,0.0) i.e in continuous feed - only BTX concentration is 
disturbed (overall BTX concentartion is not changed). Benzene is extracted up to 
25th stage (see figure 4.11[aj) from naphtha feed. BTX concentration profiles are seen 
with appropriate changes in figures 4.11(a,b) with change in feed contents of BTX. 
Dynamic holdup does not change since no feed flow rates have been changed as seen 
in figure 4.12. 

The above four logical tests (Run 2,3,4 and 5) have been performed to validate 
the model as well as computer program written for this purpose. Results obtained 
in above 4 sets support the validity of the model for the extraction column using 
the test system 1 chosen for study. Again it must be noted that, these may not 
tally with the exact values obtained from real plant as model parameters have not 
been fine tuned against real data. Similar runs have been performed on Test System 
2 (water-acetone-toluene) and found satisfactory. Results have been presented in 
Appendix D. 
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Figure 4.9: Solutes (BTX) concentration profiles (a) Continuous naphtha phase (b) 
Dispersed sulfolane phase. (Run 3) S/F increased 
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4.3 Dynamic simulation results 

Test SYSTEM 1 

Same system (Naphtha-BTX-Sulfolane) has been chosen for dynamic simulation test- 
ing and analysis. DASSL method has been used to solve set of 1440 DAE, formed 
for this system for 40 stage extraction column discussed in section 4.2. Out of 1440 
equations 480 equations are differential equations and rest are algebraic equations. 
DASSL evaluates the derivative of independent simulation variables and the set of 
1440 equations are solved by using N-R method at any instant of time; termina- 
tion criteria being same as used in steady state simulation( 4.2). Following results 
are obtained from the dynamic simulator. Four simulation runs (Run 6 to 9) have 
been performed to study dynamic behavior. Summary of these results is presented in 
Table 4.2 


Table 4.2: Dynamic simulation runs 


Run 

no. 

Total Run 
time 

Perturbation/Remarks 

6. 

60 minutes 

at 90 seconds; Naphtha feed flow rate (U in ) increased from 2.0 
moles/s to 3.0 moles/s 

7. 

90 minutes 

at 90 seconds; Sulfolane feed (l/ n ) increased from 4.5 moles/s to 
5.0 moles/s 

8. 

10 minutes 

at 50 seconds; BTX concentration in naphtha feed has been 
disturbed from (0.05,0.16,0.05) to (0.09,0.10,0.07) i.e. keeping 
overall aromatics content constant. 

9. 

2.5 hours 

at 30 seconds; Ui n from 2.0 to 3.0 moles/s 
at 600 seconds; Vi n from 4.5 to 5.0 moles/s 
at 1800 seconds; BTX concentration in naphtha feed changed 
from (0.05,0.16,0.05) to (0.09,0.10,0.07). 

at 3600 seconds; Unperturbed all the 3 perturbation given so 
far i.e. resuming feed conditions corresponding to initial state 
at time t=0 seconds. 

Note 

: In all simul 

lation runs 6 to 9, initial condition at t=0 sec has been assumed 


at steady state and corresponding to Run 1 (section 4.2). 


Simulation run 6 (for 60 minutes) (Figures 4.14 to 4.21) 


Perturbation is given in the continuous phase flow rate at 90 seconds. Continuous 
phase flow rate increased from 2.0 moles/s (0.2 kg/s) to 3.0 moles/s (0.3 kg/s) and 
data are plotted against time scale for different stages selected for analysis. Although 
each data is available for each stage at every time instant, 3-D plots can be generated 
and these are discussed later. For the purpose of analysis 2-D graphs are plotted. Fig- 
ures 4.14(a,b,c), represent the transient profile of concentrations of BTX respectively 
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in dispersed sulfolane phase. Similarly figures 4.15(a,b,c) show transient response of 
BTX concentrations in continuous naphtha phase. Concentration of BTX (solutes) 
increases with time (after 90 seconds) both in dispersed and continuous phase. Since 
dispersed phase flow rate is not perturbed dynamic holdup plotted in Figure 4.16 
registers a slight change at different stages. Step change in continuous phase flow 
rate at 90 seconds, shown in Figure 4.17 affects the continuous and dispersed phase 
flow rates exiting from the column, as seen in Figures 4.17 and 4.18. 

Simulation run 7 (for 90 minutes) (Figures 4.19,4.20,4.21) 

For this run perturbation in dispersed phase flow rate is given at 90 seconds. Flow rate 
of sulfolane phase increased from 4.5 moles/s (0.54 kg/s) to 5.0 moles/s (0.64kg/s). 
Transient response of solute (BTX) concentrations are plotted in figures 4.19(a,b,c) 
for dispersed phase and in figures 4.20 (a, b,c) for continuous phase. From all these 
response it can be seen that benzene has been extracted in less time as compared 
to toluene and o-xylene by sulfolane. This is clear in Figure 4.19(a) which shows 
benzene concentration profile has achieved steady state faster than other two solutes 
concentration profile. Figure 4.21 shows the transient response of dynamic holdup 
at four different stages, which shows significant change in its value due to change in 
sulfolane flow rate. 

Simulation run 8 (for 10 minutes) (Figures 4.22,4.23) 


For this simulation run, BTX concentration is perturbed in continuous phase feed at 
50 seconds. Naphtha feed concentration changed from (0.37,0.37,0.05,0.16,0.05,0.0) 
to (0.37,0.37,0.09,0.10,0.07,0.0), i.e. overall aromatics (BTX) concentration is kept 
constant. Effects of changing concentration in feed are plotted in Figures 4.22(a,b,c) 
for dispersed phase and in Figures 4.23(a,b,c) for continuous phase. From all these 
six plots, it can be seen that change in BTX concentration significantly affect near 
the naptha feed point (i.e. at the bottom of the extraction column) since most of the 
extraction takes place at the bottom of the column. As continuous phase moves up- 
ward, change in concentration corresponding to change in feed is not very significant; 
this confirms the early extraction of aromatics by sulfolane. Similar studies can be 
carried out for various solvents and based on their response, it can be judged which 
solvent is best suited for a particular extraction process. 


Simulation run 9 (for 2.5 hours) (Figures 4.24 to 4.34) 


In the absence of any real data, simulation package written must be checked by an- 
alyzing some logical tests. For validating the dynamic dynamic simulator, following 
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Figure 4.16: Dynamic profile of dynamic holdup (Run 6) 



Figure 4.17: Dynamic profile of forward flow (continuous naphtha phase) (Run 6) 
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strategy is devised — three (3) perturbations given in the column feed at different 
times are : 


1. At 30 seconds — Change in continuous phase flow rate from 2.0 moles/s (0.20 
kg/s) to 3.0 moles/s (0.30 kg/s). 

2. At 600 seconds — Next change in dispersed phase flow rate from 4.5 moles/s 
(0.54 kg/s) to 5.0 moles/s (0.6 kg/s). 

3. At 1800 seconds — Next change in concentration of BTX in naphtha feed from 
BTX (0.05,0.16,0.05) to (0.09,0.10,0.07), keeping other concentrations constant. 

4. At 3600 seconds — Unperturbed all three perturbation given so far. 

These three perturbations are introduced at times mentioned and after 1800 sec- 
onds (30 minutes) simulation run is further allowed for another 1800 seconds i.e total 
up to 1 hour without giving any perturbation in last half hour. After 1 hour (3600 
seconds), fourth perturbation is introduced in such a way which unperturbed all the 
previous three perturbations, thus maintaining feed conditions similar to the feed 
conditions at time t =: 0 seconds. Simulation run is further allowed for 1.5 hours i.e 
a total of 2.5 hours, in order to get steady state. The new steady state reached af- 
ter 2.5 hours matched the initial state at time t=0 seconds (Figures 4.24 to 4.34), 
also intermediate behavior is as per the expectation for any extraction process. This 
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Figure 4.19: Dynamic profile of solutes (BTX) in dispersed sulfolane phase (Run 7 ) 
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Figure 4.20: Dynamic profile of solutes (BTX) in continuous naphtha pahse (Run 7) 
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Figure 4.21: Dynamic profile of dynamic holdup (Run 7) 

confirms the validity of dynamic part of the simulation code as well as DASSL code. 
Simulation program gives data for 3-D plots i.e. value of every variable at each time 
instant for every stage. From these data transient responses have been plotted. Ev- 
ery such plot confirms that the final steady state reached after 2.5 hours is same as 
initial steady state at t=0 second. Figures 4.30(a,b,c) shows transient response for 
BTX concentration at different stages in dispersed phase and Figures 4.31(a,b,c) in 
continuous phase. Figure 4.32 plots the transient response of dynamic holdup at four 
different stages. It can be seen that perturbations get absorbed quickly by dynamic 
holdup whereas concentration profiles indicate that perturbation lasts for longer time. 
This is due to the fact that mass transfer process is slower than physical process. Sim- 
ilarilv perturbation lasts for lesser time in forward and backward flows of two phases 
inside the column. Figures 4.33 and 4.34 reflect such behavior in forward flows of 
continuous and dispersed phase flow. 

Again dynamic simulation is tested in the absence of real plant data, which may 
predict different values than shown here, but will show similar trends as described 
here. Model parameters are yet to be fine tuned for real simulation. 

Test SYSTEM 2 Similar test are performed for system 2 (water-acetone-toluene) 
and results are found similar to described above, which confirms the validity of model 
and simulation code in general (refer to Appendix D). 
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Figure 4.22: Dynamic profile of solutes (BTX) in dispersed sulfolane phase (Run 8) 
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Figure 4.25: Dynamic profile concentration in continuous naphtha phase (Run 9) 
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Figure 4.30: Dynamic profile of solutes (BTX) in dispersed sulfolane phase (Run 9) 




MOLE FRACTION MOLE FRACTION MOLE FRACTION 


4.3 Dynamic simulation results 


61 




Figure 4.31: Dynamic profile of solutes (BTX) in continuous naphtha phase (Run 9) 
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Figure 4.33: Dynamic profile of forward flow (continuous naphtha phase) (run 9) 
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Figure 4.34: Dynamic profile of forward flow (dispersed sulfolane phase) (run 9) 

4.4 Effect of axial dispersion 


One of the typical phenomena in extraction process is backmixing in flow streams 
between adjacent stages. In present model this backmixing is represented by the 
contribution of axial dispersion which causes backward flow in both continuous and 
dispersed phase flow streams as presented in section 2.2.2. Extent of backmixing dif- 
fers for different types of extraction columns and also depends upon size of equipment 
and operating conditions. The backmixing effect would be high in agitated columns 
such as RDC, pulsed sieve tray columns etc. , because of high turbulence. However 
such effects are not significant in unagitated columns such as spray columns and un- 
agitated sieve tray column. Since we are considering an unagitated sieve tray column, 
backmixing is assumed only in continuous phase. Backmixing in dispersed phase is 
neglected completely by the fact that separated compartments of the column by sieve 
tray may not allow mixing of dispersed phase fluid in adjacent two compartments. 
Backmixing effect is represented by the axial dispersion coefficients (D c , D d ) in model 
presented here 3 . Correlations of these dispersion coefficients for few types of columns 
are available in literature (Baird [4]) however these correlations are not well tested 
and reliable enough to be incorporate in this study. 


Simulation run 10 (for steady state) (Figures 4.35 to 4.39 

3 Axial dispersion with in stage fluid may be high enough, which supports complete mixing 
assumption 
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Figure 4.35: Effect of axial dispersion on solute (BTX) profile in dispersed sulfolane 
phase (A) D c = 3cm 2 /s (B) D c = 35 cm 2 /s 

In order to depict axial dispersion affect on extraction process, fixed values of disper- 
sion coefficients D c = 3.0cm 2 /s and D d = 0.0 are taken for one simulation run and 
D c = 35cm 2 /s and D d = 0.0 for another simulation. Results obtained from both 
simulations are imposed together and plotted in Figures 4.35 to 4.39. Higher axial 
dispersion 4 , should cause loss in extraction efficiency. This fact is supported by the 
results obtained in Figures 4.35 and 4.36. Figure 4.36 shows less extraction from con- 
tinuous phase, especially of toluene. Figure 4.37 shows the dynamic profile holdups 
for these two runs. Figure 4.38 shows effect of change in D c value on continuous 
phase forward flow rate sudden change in flow rate is due to higher value of axial 
dispersion coefficients chosed for stages but not at the ends. Similar change is seen 
from Figure 4.39 for dispersed phase. These trends are matching with the results 
published by Zimmermann [7] for pulsed sieve tray column. 



4 means higher backmixing 
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Figure 4.37: Effect of axial dispersion on dynamic holdup profile 
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Figure 4.38: Effect of axial dispersion on forward flow (continuous phase) 



Figure 4.39: Effect of axial dispersion, forward flow (dispersed phase) 
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4.5 Model parameters 


Before going to drop breakage and coalescence effects on extraction process, it is 
necessary to discuss about model parameters which have to be fine tuned by opti- 
mization against real operation data of extraction column. Model used in this study 
is developed incorporating maximum theoretical considerations and minimum of free 
parameters. But due to lack of perfect theoretical knowledge about axial disper- 
sion, multicomponent mass transfer coefficients and drop breakage and coalescence, 
empirical correlations have been used to fill this gap. 

For axial dispersion correlation have been developed by several authors ( Baird [4, 
pp 323]), but they are very specific to column size and type of a column. Axial 
dispersion coefficients D c and D d may therefore be taken as model parameters which 
have to be tuned for a particular extraction column before applying present model 
for that column. 

Correlations used for calculating dispersed phase velocity u d and relative veloc- 
ity u; (see equation B.4 section B.1.2) have empirical constants. Similarly correla- 
tions used for drop breakage and coalescence such as breakage frequency, collision 
frequency and coalescence efficiency (see eqn. B.5 to B.9 ) also contain empirical 
constants. These should be determined by experimentation and optimization there- 
after. Several attempts have been made in literature to evaluate the values of these 
empirical constants by carrying out independent studies 5 . The effect of mass transfer 
may not be significant on these empirical constants; by assuming this, same empiri- 
cal constants have been used for the mass transfer model in this thesis. However if 
some of these empirical constants show significant change in their value when mass 
transfer effects are incorporated then these empirical constants should be included in 
the list of model parameters and must be determined simultaneously by conducting 
hydrodynamics cum mass transfer experimentation and optimization thereafter. 

In this thesis effort has been made to evaluate mass transfer coefficients rigorously 
using multicomponent mass transfer theory as developed by Taylor [21]. In some pre- 
vious works, such as Spencer [3] and Ricker [2], binary mass transfer coefficients have 
been considered as model parameters and evaluated by optimization for a particular 
column and extraction process. In present model if binary mass transfer coefficients 
are not supplied then they are calculated from empirical correlations (equations B.10 
to B.13) for particular column and then interactive mass transfer coefficients are 
evaluated. The use of such correlation may introduce some errors in final result. 

Thus only D c and D d can be taken, in the first selection as model parameters. 
In case of dissatisfatory results other parameters the empirical constants mentioned 
above may be included in the list of model parameters. 


r> generally without considering mass transfer 
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4.6 Effect of drop breakage and coalescence 


As mentioned in the previous section correlations used for breakage frequency g(dj), 
collision frequency h(di,dj) and coalescence efficiency A (d{,dj) have been tested for 
agitated columns in literature. Empirical constants used in these correlations have 
been optimized for agitated columns. Values of these constants used for our work 
are mentioned in Appendix B with appropriate references to the type of column they 
were determined for. Since same values have been used in our work for unagitated 
sieve tray column the results mentioned in this section may not be highly reliable. 
For the sake of completeness simulation runs have been performed by taking drop 
breakage and coalescence into account for unagitated sieve tray column. 

In order to study only drop breakage and coalescence effects only 4 stage with 
8 drop classes have been considered for system 1. instead of 40 stages considered 
in earlier sections. For unagitated sieve tray column drop spectrum (minimum to 
maximum diameter of drops) may not be very wide. Following 8 drop classes have 
selected (refer to section 2.8.2) by the simulator after discretization of drop spectrum. 


Drop class 

Diameters (mm) 

(l) 

1.43 

(2) 

1.82 

(3) 

2.08 

(4) 

2.29 

(5) 

2.46 

(0 

2.62 

(7) 

2.76 

(8) 

2.88 


Four simulation runs are carried out as 

(a) With drop breakage and coalescence. 

(b) With only drop breakage. 

(c) With only drop coalescence. 

(d) Without drop breakage and coalescence. 

Holdup profiles for all the eight drop classes are indicated for these four runs 
in Figures 4.41(a,b,c,d). Total holdup for these four runs has been compared in Fig- 
ure 4.42. Corresponding Sauter mean diameter of drop swarm in a stage are plotted in 
Figure 4.43. No significant effect of drop breakup and coalescence is noticed in these 
plots, since drop spectrum does not change much due to weak drop breakage and 
coalescence in unagitated sieve tray columns. Sauter mean diameter show decreasing 
trend for breakage but increasing trend for drop coalescence as seen in Figure 4.43. 
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Figure 4.40: Key for Figures 4.41 
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4.41: (a) With drop breakup and coalescence, (b) only drop breakage, (c) only 
>aleseence. (d) without drop breakup and coalescence. 
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Chapter 5 


Conclusions and Suggestions 


Simulation model selected for present thesis incorporates non-equlibrium approach 
and drop breakage and coalescence model to describe dispersed phase behavior. The 
model has been extended for dynamic simulation in this thesis. In general results, 
predicted by the model are found satisfactory and as per the expectations. Model 
has been tested for two extraction systems — Aromatic extraction by sulfolane and 
acetone extraction by toluene. 

Axial dispersion coefficients for continuous and dispersed phase are assumed as 
model parameters and given some suitable fixed value during simulation. 

For distribution coefficient calculation UNIFAC method has been used. For LLE 
UNIFAC data at 25° (7 are available in literature for all of the components which have 
been used in this work. However extraction column described for System 1 works at 
an average temperature of 95°C 1 . This temperature difference may cause error in 
final results which could be corrected by taking the UNIFAC parameters at 95°C r . 
Modified UNIFAC method may be used which includes temperature effects also but 
data published in literature are limited for this method. 

One of the important aspect of present model is its non-equlibrium approach to 
deal with interphase mass transfer and other is to incorporate drop breakage and 
coalescence for dispersed phase by using drop population approach. Effects such as 
osmotic diffusion and reverse diffusion are accounted for in the rigorous analysis used 
in present study to calculate interactive mass transfer coefficients. The validity of 
these results will definitely increase if experimental values of binary mass transfer 
coefficients are available; otherwise a suitable correlation can be used to predict de- 
sired estimates as in this thesis. A detailed study on interaction of species during 
interphase mass transfer is now possible by using such model presented here. 

For unagitated column studied in this thesis drop breakage and coalescence create 
no significant changes in process variables. Further, the empirical constants used are 

1 Mukhopadhyay [28] has given UNIFAC data at 100 ° C but not found satisfactory with this work 
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those reported in literature for agitated columns. Especially coalescence efficiency is a 
very critical parameter which is not well understood. During steady state simulation 
its value has been evaluated in the range of 5-15% for correlation used for our work. 
Some other correlations Tsouris [17] predict even a very low value (less than 1%) 
of coalescence efficiency for drop diameters used in System 1. One has to take a 
judicious decision to select proper correlation. 

The motivation behind present work is to develop a suitable model for liquid- 
liquid extraction which can be used for detailed study such as required for the design 
as well as for the control of extraction process. One of the major advantage of our 
formulation is that it can be applied for any type of countercurrent extraction column 
once sufficient information on design parameters and correlations is available for that 
column. For drop breakage and coalescence, models used ealier without considering 
mass transfer effects have been modified for mass transfer effects and included with 
the model presented in this work. Such an inclusion now allows study and control of 
physical parameters such as holdup in conjunction with mass transfer. 

One of the crucial decisions in simulation is to select a robust algorithm for solving 
set of equations. Newton-Raphson method for solving set of nonlinear equations has 
been found robust in this study. Significant amount of time and computer memory 
saving has been obtained using Block-Thomas algorithm to solve the set of linear 
equations 2 in which block tridiagonal sparse matrix is formed. Smooth convergence 
of equation solver was also our major concern which has been achieved successfully 
by using N-R method. Practically, evaluation of Jacobian in every iteration takes too 
much time, especially evaluation of numerical derivatives. Finite difference approx- 
imation of numerical derivatives greatly increases the cost of solution because more 
physical properties evaluations are required; further, neglect of these derivatives may 
increase number of iterations or even cause failure. For dynamic simulation differen- 
tial algebraic equation (DAE) set is formed. For solving DAE , DASSL which is one 
of the robust methods has been selected in this study. DASSL has certain advantages 
over other methods such as Runge-Kutte (RK) and Gears method. Generalized RK 
method requires large computer memory and longer time for solving DAE problem. 
DASSL is a variable order and variable step size method; whereas in Gear’s method 
using backward difference formula (BDF) formula the order is fixed. The selection 
of suitable order in DASSL during simulation allows large time step thus reducing 
time for dynamic simulation. On the other hand Gear’s method having a fixed or- 
der resorts to small step size in a similar situation which slows down the simulation. 
Another advantage of selecting DASSL in this work is that it is simple and easy to un- 
derstand; it has been found logistically convenient to incorporate DASSL with steady 
state simulation computer code. However for RK method or Gear’s algorithm major 
part of the dynamic simulation problem had to be reformulated. DASSL has been 
tested on variety of simulation problem such as simulation of mechanical, electrical 
and aerodynamic systems satisfactorily. 


2 formed after linearizing, using N-R method 
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Present model can be studied in depth of rigorously by including several drop 
classes, mass transfer coefficients, drop breakage and coalescence affect and UNEFAC 
method. On the other hand it can be studied by neglecting or simplifying any of 
these details to reduce time of computation with out losing much on accuracy. 

SUGGESTIONS FOR FUTURE WORK 
Modeling of any chemical process is a creative exercise and success of modeling is 
complete only when the model is tested against real operational data. Feedback from 
real process may further suggest some refinement in modeling aspects which can't be 
envisaged just by carrying out theoretical study. It is therefore recommended that 
present model must be tested for a variety of extraction systems which should help 
to remove any bug and make the simulator versatile for study of any new extraction 
process. Other suggestions to make add on to the effort in this thesis, for a better 
success in future are : 


1. Model parameters such as axial dispersion coefficients have to be tuned by 
optimization against real process data. Tuning of empirical constants used in 
correlations may be done for unagitated columns by carrying out independent 
experimentation as done in literature for agitated columns or they can be tuned 
with present model. However latter option increases the dimensionality of the 
optimization problem. 

2. Model should be tested for agitated columns. For these columns abundance of 
data are already available in literature. Computer program written as a part of 
this thesis can be modified for such study by replacing correlations for agitated 
column. 

3. In order to reduce CPU time 3 for solving set of non-linear equations, hybrid 
techniques such as Quasi-Newton method and combination of Newton and 
Quasi-Newton suggested by Lucia [44, 45] can be used. These methods up- 
date Jacobian without recalculation of each element - especially the Jacobian 
part which contains numerical derivatives is updated by the algorithm proposed 
by Lucia. Instability however may creep in by using such techniques. 

4. In our work dynamic simulation is carried out from one steady state to another 
steady state i.e. DASSL is applied at initial steady state and then perturbations 
are introduced after some time steps. It would be useful to extended this for 
startup and shutdown of extraction column also. In our work time derivatives 
at time t=0 are taken zero (y'{ 0) = 0). DASSL can also be applied during 
transient state for which time derivatives have to be supplied or evaluated (see 
Brenan [37] for further details). 

3 One can use “prof’ command in UNIX to get time consumed by each subroutine during 
simulation. 
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5. Though isothermal operation is assumed in our work, sometimes temperature 
profile is appled along the column height. In such case same model can be used 
by assuming any stage at constant temperature, however temperature difference 
may exist from one stage to next. Computer program has been made in a 
modular style keeping this factor into account. Coding of physical properties 
relations incorporating temperature dependence in appropriate subroutine will 
facilitate the non-isothermal simulation. 

6. The detailed non-equilibrium dynamic model using population balance ap- 
proach is not suitable as a process control model. Reduction of model by 
clubbing several real stages into one stage, using fixed distribution ratio (K- 
value) and fixed mass transfer coefficients for a range of operating conditions 
and avoidance of drop breakage and coalescence can be resorted to carry out 
modeling for control. 

7. Major improvement in computation time can be achieved if some equations 
as well as independent variables are eliminated by substitution. For example 
equilibrium relations EQijk and transport relations MBIDij k , MBIC^ can be 
clubbed together by substituting for either x\ ik or y- jk from equilibrium relation 
into transport relations. This will reduce size of the simulation problem and 
memory requirement also. Such change will require reformulation of Jacobian 
and consequently some part of computer code also. 

8. Drop breakage and coalescence study can be extended further and some useful 
results may be obtained for different types of column. Presently this field is 
under intense research in the field of liquid-liquid dispersion. Most of the re- 
lated work published in literature does not include mass transfer effect. Model 
presented here has been modified specifically for studying drop breakage and co- 
alescence under mass transfer conditions. A laboratory scale extraction column 
is recommended for such study. 
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Appendix A 

DASSL : Differential Algebraic 
System SoLver 


This section describes briefly about DASSL formulation and its implementation in 
this study. Details of this method and some useful definitions commonly used for 
solving Differential Algebraic Equation (DAE) can be seen from Brenan, Campbell & 
Petzold [37, 1989]. DASSL is a very robust tool for solving DAE of index 1 0 and 1. 
It also solves the equation set of semi-implicit type of DAE (as is in this study) very 
well. DASSL requires the equation set to be solved in the form of 

E ^tn+l> t/n+l: Un+l) = 0 (A.l) 

with 

y(t 0 ) = y 0 

y'{to) = y'o 

where h n+ i = t n+ i — t n step size at current step (n— 1). The non-linear set of equa- 
tions A.l is then solved by using multivariable Xewton-Raphson method. In this 
study 2/' (to) is evaluated by assuming steady state at time t=0, i.e. y'(t 0 ) = y' 0 = 0 is 
taken. This is true for one steady state to next steady state but may not be true for 
startup/shutdown problem. DASSL approximates the derivative y' n+1 using the k th 
order backward differentiation formula (BDF), where k ranges from 1 to 5 (suggested 
by Brenan [37, pp 118]). DASSL is a variable order and variable step size method 
which makes it most robust and powerful to solve even steep differential equations. 
Detailed derivation can be seen from Brenan [37]. Here a brief formulation of DASSL 
is given : 


1 Index of DAE defined as number of times differentiation required to convert DAE into ODE set. 
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A.l Intermediate parameters used in DASSL for- 
mulation 


Some intermediate parameters used for DASSL are : 


i>i(n+ 1) 

cti(n + 1) 

Pi(n + 1) 

Pi(n + 1 ) 

M n ) 

M n ) 

cri(n + 1 ) 

<Ji(n + 1 ) 

7 i(n + 1) 
7 i(n + 1 ) 

Qtg 

a°(n + 1 ) 


hn + 1 “I” h n + * • • + h n +2—i — ^n+1 Ai+1— i ? 2 1 

= h n+ i/tpi(n + 1) 

= 1 

_ ^i(n + VyMn + 1) • • •‘0j- 1 (n + 1) . 

^i(w)^(n) ••••&_! (to) ’ % 

= Vn 

= .[y n y n _ i+1 ] , »>1 

= /3i(n + l)0*(n) , z > 1 

= 1 

_ ^n+l (i ~ 1 ) ! j > ]_ 

•0i(n + l)'^ 2 (« + l)---'0i(n + l) 

= 0 

= li-\{ n + 1) + + l)/h n+ i , i > 1 



k 

= ~X>j> + i) 

i=i 


(A-2) 

(A.3) 

(A-4) 

(A.5) 

(A.6) 

(A.7) 

(A.8) 

(A.9) 

(A.10) 

(A.11) 

(A.12) 

(A.13) 

(A-14) 


The divided difference 2 [y n ,y n _ lt . ■ - ,y n _ i+1 ] formula can be seen from any standard 
numerical method text. 


A. 2 DASSL formula for time derivatives y' n+l 


DASSL formula for time derivative of variable is given as : 

y'n+l =Yn+l - r^”(yn+l - Yn+l) 


where, 


and 


J-n+l 
fc+l 

^n+l = E©i( n ) 

i=l 


fc+1 

y'n+l = Y,°/i( n+ 1 ) ( ^*( n ) 

i= 1 


(A.15) 

(A.16) 

(A.17) 


“eg: [»l,So] = *?+§ 
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A. 3 Step size selection 

In this study following strategy is adopted to decide the correct step size. At any 
step, for correct step size and order, solution of DAE must satisfy following error 
condition. 

ERR = M\\y n+ i — y° +1 || < 1.0, (A.18) 

where 

M= max (a fc+1 (n + 1) , | a fc+1 (n + 1 ) + a s - a°(n + 1) |) (A.19) 

• At step=l, using ho = 10 -3 x (TOUT — 0) as previous step size for time step 

1. Here TOUT is the time for which simulator has to be run. 

• Take successful step size of previous step, raise the step size in next step, if it 
can be doubled 3 , h n+ i=2h n , solve DAE and check error condition A.18 

• If raised step size is not accepted then use previous successful step size, h n+ i=h n . 

• If above two steps failed then step size is to be decreased, follow the sequence if 
this decrease is called for. 

1. Step size is decreased at least by a factor r=0.9 i.e. (/i n+1 =rh n ) and at 
most r=0.5 and see where it could be accepted; if not, call it failure no. 1. 

2. If failure 1 occur, then take r again from 0.9 to 0.25 and multiply by the 
step evolved after failure no. 1; keep on checking error condition after 
every new step selected; if not satisfied, call this as failure no. 2. 

3. After failure no. 2 the step size is reduced by one quarter (1/4), based on 
the philosophy that such a high step is not trusted. 

• Keep on counting the numbers of failure in above operation. After a large step 
size decrease (i.e. too many failure), the BDF formula is most accurate at order 
one. Therefore reduce current order to 1 and keep on trying again; after each 
failure reduce step size by 1/4 till [ t + h « t ], where it is declared that DASSL 
FAILED. 


A. 4 Order selection for next step 


Initially order is kept at 1, and step is increased as described in section A.3. After 
that the order is raised or lowered depending upon whether error estimates of Taylor 
series - TERKM2,TERKM1,TERK and TERKP1 form an increasing or decreasing 

3 Sometimes this doubling may give severe instability problem in equation solver so vary this 
factor from 1.5 to 2.0, try to keep as high as (upto 2.0) possible 
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sequence. If an increasing sequence is formed then order is reduced and if a decreasing 
sequence is formed than order is increased. The error estimates are given below. 


TERKM2 = ||(fc - l)er A ._ 1 (7i + l)<f> k (n + 1)|| (A.20) 

TERKMl = pcr fc (n + l)<£ fe+1 (n+ 1)|| (A.21) 

TERK = \\(k + l)a k+1 (n + l)<f> k+2 (n + l)\\ (A.22) 

TERKP1 = \\(k + 2)<7k+2{ n + l)^fc+3( w + 1)|| (A.23) 


The DASSL technique is still under refinement, Hence in case of failure some 
suitable strategy has to be adopted based on intuition, which may be specific to the 
problem. If process behavior or its estimate is available a priori some time step saving 
steps could be adopted. In this study continuing with DASSL, if any perturbation 
is given to the process then simulator (computer program) integrates DAE first up 
to the time of perturbation. After that very small step size is taken (irrespective of 
what DASSL suggest) to introduce perturbation and finally DASSL is again allowed 
to implement step size and order as described in sections A.3 and A.4. 



Appendix B 
Correlations Used 


The model described in this study is general and could be used for agitated as well 
as unagitated multicomponent countercurrent extraction column . However, in this 
study it is tested for Unagitated Sieve Tray extraction column. For other columns 
same computer code can be used with modifications in the correlations used. Follow- 
ing sections present the correlations used to test the model for unagitated sieve tray 
column. 


B.l Design parameters 

B.l.l Drop diameters 

• Maximum stable drop diameter d max 
taken from (Baird [4],pp 129, eqn. 17) 

dm “ = c 'i^) (e) "°' 4 (S- 1 ) 

where C\ is a constant depending upon system properties and type of the col- 
umn. For RDC value of C\ = 0.72 is reported. For this study its value is taken 
as 0.8 considering the fact that in an unagitated column larger drop diameter 
may exist. 

e is power dissipated per unit mass of fluid phase. Relationships of e for some 
contactors are given in (Baird [4], pp 130, Table 1). For our work e relation 
for spray column, which is assumed to be applicable for each compartment of 
unagitated sieve tray, is taken : * 

e = u i9 ^- (B.2) 

Pc 

Ud is a superficial velocity of dispersed phase. 
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• Mean diameter d mean or d 32 or d vs 

Several reliable correlations are available for mean drop diameters possible in 
any contactor. Various authors call it by various names such as dz 2 , vediyan 
mean drop correlation d vs etc. The correlation proposed by Vediyan (Laddha [1], 
pp 308, eqn. 6.53) has been used for unagitated sieve tray column : 


d v 


{~f/Apg) 


1.592 



— 0.0665 


(B.3) 


where U n and d n are jet velocity and nozzle (hole) diameter respectively. 


B.1.2 Characteristic velocity u* 


Various author have published several correlations. The correlation proposed by 
Vediyan (Laddha [1], pp 313, eqn. 6.54) is used here. 

- 0.0818 

(B.4) 

Other correlations are also available from Baird [4], Khani [8]. 


ul 


(ta pg/piy 


= 1.088 


'J£ 

2 g a 


B.1.3 Static holdup, Jet velocity, Jet diameter . . . 

These are used in design calculations. All these calculations are referred from (Lad- 
dha [1], Chapter 11- Perforated Plate Towers). These calculations are specific to 
perforated towers only and therefore not valid for other types of column. 


B.2 Drop breakage and Coalescence 

B.2.1 Breakage frequency g(dj) 

Mechanistic model (which includes the physical properties of the system ,the geometry 
and energy provided to the dispersion) was first proposed by Coulaloglou 10] in 19/7 
for agitated columns. Relations cited by (Tsouris and Tavlarides [17], pp 396, eqn. 
4) are used here : 

(B.5) 

wh^re k \ , k 2 are empirical constants, evaluated by experimentation and optimization 
for a particular column. Hsia [13] evaluated these parameters for dispersion in stirred 


g(dj) = ki 


1 

63 


(1 + <t>)dj 


—exp 


-h 


7(1 + 4 >f 

! i 5 

Pd&d] J 
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tank. In the absence of any information of these data for unagitated sieve tray 
column, same values of these parameters have been used for this study also. It 
has been assumed that order of these parameters may be same for other columns. 
However experimentation is necessary in order to validate results for drop breakage 
and coalescence for other extraction columns. Typical values are k x = 0.00481 and 
k 2 = 0.0551 taken from Hsia [13]. 


B.2. 2 Daughter drop distribution 


After breakage of any drop the produced drops are called as daughter drop. It is 
assumed that father drop after breaking gives certain spectrum of daughter drops. 
The probability density of this spectrum is called as daughter drop density function 
represented by (3(di,dj). It represents the probability density of daughter drop of 
diameter dj due to breakage of father drop of diameter d z . For convenience sometimes 
it is represented as or /3(vi,Vj) also. The distribution used in PhD thesis of 

(Tsouris [16], PhD thesis, pp42, eqn. 4.20) is used here : 


P(v, v') 


a 


(” - 


SL, {v-’ifdv 


4 “ 


a 


V - Vr, 


where a is given by 


a = exp 


-c 2 - 


d dnn 


dn 


(B.6) 


(B.7) 


Above c 2 is a constant the optimum value of which has been found to be 0.75. Critical 
or minimum drop diameter d m i n is taken as 0.0 for present study. The value of 
parameter a lies between 0 and 1, which are its extreme values. 


B.2. 3 Collision frequency 

Theoretical formulation can be seen from Kennard [39, 1938.pp 103-114], in where 
collision frequency formula has been developed for gas molecules. The formulation 
used here is developed by (Tsouris [17], pp 401, eqn. 50). 

h(d h dj) = ^{di + djfiuf + u^riirij (B.8) 

where mean velocity is given by it? = 1.07eidt. e here used for spray column as in 
equation B.2. 

B.2. 4 Coalescence efficiency 

This is not well understood in the literature published so far. Some definition favour 
higher efficiency for smaller drops and some for larger drops . For present, work model 



B.3 Binary mass transfer coefficient 


86 


proposed in PhD thesis by Tsouris [16] is used. 


A (d{, dj) = exp 


c Uh 

. (dj + dj)$ _ 


(B.9) 


where optimum value of constant C3 has been found to be 8.0 x 10 2 for a column 
used in Tsouris’ thesis. In the absence of any data same value is used here. Typically 
coalescence efficiency ranges between 5-15 % using above formula. 


B.3 Binary mass transfer coefficient 


Although it is preferd to use experimental values of binary mass transfer coefficient, in 
the absence of such data, correlations have to be used for calculations. These correla- 
tions are geometry dependent and are the basic inputs for the theoretical rigorous cal- 
culation of interactive mass transfer coefficient matrix used in MBIDijk and MBICijk 
equations sets. For present study, correlations Tables published in (Baird [4], pp 324- 
326, Tables 1,2) are referred for binary mass transfer coefficients. The simulator has 
been tested using the following two correlations for continuous and dispersed phases 
respectively. 


kj 


d 

0.0432- 

1/ 


Pd A 
.Md) av 


\dg j 



Pd \ 

Vpddr/J 


0.601 


(B.10) 



(B.ll) 


These are used for mass transfer coefficients prediction at the time of drop detachment 
from sieve tray. Rise time t is given by 


n 0 (7r/6)d 3 

Qd 


(B.12) 


Further {jfc) an d (f^) are equvivalent to molar densities c d and c c respectively. 
For diffusivity calculation famous empirical correlation by Wilke-Chang (1955) is used 
see (Taylor [21], pp 91, eqn 4.2.18,4.1.8) or (Reid [42], pp 549, eqn 11-45) and is given 
below. 


V im = 7.4 x IQ’ 8 


/ (0 m M m ) l / 2 r \ 
V Pl^Vf- 6 ) 


where 

V im = mutual diffusivity of solute i into m (cm? / sec) 

Mm — molecular weight. 

T = temperature K 

— viscosity ( centipoise) 

Vi = molal volume at normal boiling point. (cm?/g mole) 
(j) m = association parameter. 


(B.13) 



Appendix C 

Computer program and data 


A stand alone computer code written in “C” language for simulating multicomponent, 
countercurrent staged extraction process. The program simulates for both steady 
state and dynamic simulation. Complete simulation package contains 19 computer 
files with 12,000 lines of “C” code. This program has been tested on UNIX machines 
like DEC 2000 aXP and HP9000/850. Following is the list of the files and brief 
description about each file. 


Table C.l: Computer files in simulation package developed. 


S.no 

Computer file 

Brief description 

1 

simuldefs.h 

Definition file; Declares and define all functions and def- 
initions used in whole simulation code. 

2 

main.c 

Startup file, contains main program 

3 

simulations 

Selects steady state or dynamic simulation options and 
sets appropriate flags. 

4 

step.c 

Comprises main simulation loop. It calls various other 
functions like initial guess settings, conditioning input, 
N-R method step size optimization etc. 

5 

DASSL.c* 

Main body of DASSL code and formulation of dy- 
namic simulation. It also calls various functions from 
simulfun03.c 

6 

changes* 

It reads perturbations from file supplied by user which 
stores pertutbation during simulation. This file also 
made some checks before applying any perturbation. 
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7 

block.c 

The functions stored in this file are called by file 
step.c frequently for calculation of A,B abd C blocks of 
Jacobian. 

8 

vector.c 

Similar as block.c , but call for function vector (residuals) 

D vector. 

9 

a_minorx 

Calculates 30 different minors A(0,0) to A(4,5), for A 
block. 

10 

b_minor.c 

Calculates 30 different minors B(0,0) to B(4,5), for B 
block. 

11 

cjminor.c 

Calculates 30 different minors C(0,0) to C(4,5), for C 
block. 

12 

d-vector.c 

Calculates 5 different subvectors D(0) to D(4), for D 
vector. 

13 

simulfunOl.c 

Stored various functions which are used for calculating 
expressions such as average mol. wt., axial dispersion 
coefficients, molar density etc. 

14 

simulfun02.c 

Extension of simulfunOl.c, but stored most of the corre- 
lations also. 

15 

simulfun03.c 

Similar to simulfun02.c, extensively used for dynamic 
simulation only. 

16 

break_coal.c 

This stores functions for drop breakage and coalescence 
calculations. 

17 

thermo_model.c 

Stores about UNIFAC method. 

18 

mtc.c 

For mass transfer calculations. 

19 

simutils.c 

General utility function file, for programming point of 
view. Functions stored here are general in nature and 
may be used with any other “C” program. 


[*] Are the files, used only for Dynamic simulation, never by steady state simula- 
tion. 


A computer flowsheet of whole simulation package is shown in Figure C.l both 
for steady state and dynamic simulation. 
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Figure C.l: Flow sheet of computer program for simulating extraction column 


TFINAL * Total run time for dynamic simulation. 
TO a initia time (0 seconds) 
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Following are the sample data files for simulation run no. 1, described in chapter 4. 


Table C.2: Column 1 specification 


Line no. 

Data in computer file 

Description 

1. 

0.300000 

Column diameter (m) 

2. 

0.100000 • 

Tray spacing (m) 

3. 

40 

No. of real stages 

4. 

0.003175 

Hole diameter (m) 


0.010000 t 

Pitch and type of pitch (triangular) 

6. 

384 

No. of holes/plate 


Table C.3: System 1 specification 


Line no. 

Data in computer file 

Description 

1. 

6 

No. of components 

2. 

100.20 98.18 78.11 92.13 106.16 120.00 

Mol. wt. of each 

component* 

3. 

1 

No. of drop class 

4. 

0.007500 

Interfacial tension 

5. 

0.000320 0.003000 

Viscosity of continuous & 
dispersed phase 

6. 

750.0 1260.0 

Density of continuous & dis- 
persed phase 

7. 

0.000209 0.000240 0.000260 0.000290 

0.000345 0.003000 

Viscosity of each component 

8. 

0.000162 0.000140 0.000097 0.000119 

0.000140 0.000120 

Molar volume of each 
component 

9. 

1.000000 1.000000 1.000000 1.000000 
1.000000 1.000000 

Association parameter of 
each component 


DATA ARE IN SI UNITS 


(*) n-Hep- 

tane(l), cyl-Hexane(2), Benzene(3), Toluene(4), o-Xylene(5) and Sul- 


folane(6). 
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Table C.4: Operating condition 


Line no. 

Data in computer file 

Description 

1 . 

500000.0 

Operating pressure ( N/m 2 ) 

2. 

368.15 

Operating temperature (K) 

3. 

2.0 

Molar feed flow, continuous 
phase (moles/s) 

4. 

0.37 0.37 0.05 0.16 0.05 0.00 

Composition (mole fraction) 
feed Naphtha 

5. 

4.5 

Molar feed flow, dispersed 
phase (moles/s) 

6. 

0.00 0.00 0.00 0.00 0.00 1.00 

Composition (mole fraction) 
feed Sulfolane 


Next two data lines will be entered for each 
stage side streams, if they are not present 
enter 0.0 for data. First these lines are en- 
tered for continuous phase (for all 40 stage) 
and then same exercise repeated for dis- 
persed phase side streams. 


7. 

0.00 

Side feed to stage 1 (flow 
rate) 

8 

0.00 0.00 0.00 0.00 0.00 0.00 

its composition (mole 

fraction) 

9. 

10. 

•• 


11. 

so on .. 



Table C, list parameters used for UNIFAC model, which are supplied as a separate 
file. For system 1, 6 UNIFAC groups (Magnussen [24]) are identified as : 

— CH 3 , — CH 2 , — CH, — ACH, — ACC Hz and — TMS (Separate group for 
Sulfolane). 


TABLE C . 5 : UNIFAC parameters 


LINE(l) ==>(no-group) Total Number Of Groups 

LINEC2) ==>(Rk) Volume parameters (Rk) for groups (In Particular order ) 

LINEC3) ==>(Qk) Surface Area parameters (Qk) for groups (In same order) 

LINEC4-9) ==>(nu(k, i)) No. Of Groups Of Type "k" in Molecule "i". 

Enter In one LINE for One specie (molecule) , as 


nu(l,l) nu(2,l) 


nu(no-group,l) 
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Line no. 

Data 

1 . 

6 

2. 

0.9011 0.6744 0.4469 0.5313 1.2663 4.0358 

3. 

0.8480 0.5400 0.2280 0.4000 0.9680 3.2000 

4. 

250000 

5. 

151000 

6. 

000600 

7. 

000510 

8. 

000420 

9. 

000001 

10. 

0.0000 0.0000 0.0000 -114.8 -115.7 561.40 

11. 

0.0000 0.0000 0.0000 -114.8 -115.7 561.40 

12. 

0.0000 0.0000 0.0000 -114.8 -115.7 561.40 

13. 

156.50 156.50 156.50 0.0000 167.00 21.970 

14. 

104.40 104.40 104.40 -146.8 0.0000 238.00 

15. 

67.840 67.840 67.840 59.160 26.590 0.0000 


(This data set are taken from Magnussen [24]. 


nu(l,no-spicie) mi (no-group ,no-spicie) 

LINE(10-15) ==>(a(m,n)) Interaction parameters,- "m" to "n" . as 

a(l,l) a(l,2) a(l, no-group) 


a (no-group, 1) a (no-group, no-group) 


All these data files are also supplied for dynamic simulation also. Apart from these 
a perturbation file is also supplied to computer program in case of dynamic simulation, 
which stores perturbation information to be applied during dynamic simulation. This 
sample computer file for simulation run 6 (see section 4.3)is reported here as it is. 


Perturbation file (Simulation run 6) for dynamic simulation 


#Read These INSTRUCTIONS CAREFULLY !'.! 

#[1] Try Not to Edit Any Text or Insert any Text line in this file, Use . 

# [ 2 ] if required. Use # character in first column to ignore that line 

# [3] Enter Data Change Only (in place of * ),If Particular variable is not 

# Changing write * (star) in place of that 

#[4] DO NOT REMOVE ANY LINE OR CHANGE DATA SEQUENCE 

#Ex ample Entry : 
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# 

# Change No.= 5 , At Time = 45.0 seconds 

# Uin = 0.95 

# Vin = * 

# 

#Indicates at perturbation no 5, time 45.0 seconds ,Uin changes from prev. 
#to 0.95, Vin remain same 

# Only blanks and TABS axe allowed between = and data value or * 
#########################################################################* 

# 


Total No. Of Changes (perturbation) proposed = 1 

Change No.= 1 , At Time = 90.000000 seconds 

Uin = 3.0 

yin[l] = * 

yin [2] = * 

yin [3] = * 

yin [4] = * 

yin [5] = * 

yin [6] = * 

Vin = * 

xin[l] = * 

xin [2] = * 

xin [3] = * 

xin [4] = * 

xin [5] = * 

xin [6] = * 

Fc[l] = * 

yf [1] [1] = * 

yf [2] [1] = * 

yf [3] [13 = * 

yf [4] [1] = * 

yf [5] [1] = * 

yf [6] [1] = * 

Fd_in[l] = * 

xf[l][l] = * 

xf [2] [1] = * 

xf [3] [1] = * 

xf [4] [1] = * 

xf [5] [1] = * 

xf [6] [1] = * 
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Fc [40] = * 

yf [1] [40] = * 

yf [2] [40] = * 

yf [3] [40] = * 

yf [4] [40] = * 

yf [5] [40] = * 

yf [6] [40] = * 

Fd_in[40] = * 

xf [1] [40] = * 

xf [2] [40] = * 

xf [3] [40] = * 

xf [4] [40] = * 

xf [5] [40] = * 

xf [6] [40] - * 



Appendix D 


Results for System 2 


This Appendix list the results obtained for system 2 i.e. (Water (1), Acetone(2), 
Toluene(3). 



Appendix D 


Results for System 2 


This Appendix gives the results obtained for system 2 i.e. (Water(l), Acetone(2), 
Tolueno(3). This system is recommended for the lab experiments in which Toluene 
as a solvent (dispersed phase) enters from bottom of a 4 stage unagitated sieve tray 
column, which extracts Acetone from Water. Continuous phase feed enters from top 
of the column. Specifications are in Table D.l. 
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Table D.l: Specifications for for Test System 2 



200000.00 


298.15 

0.65 

0.90 0.10 0.00 
0.32 

0.000 0.000 1.000 
0.00 

0.00 0.00 0.00 


0.00 

0.00 0.00 0.00 
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Table d7i Contd.. 

UNIFAC data (as in Table 

"05)~ 

5 

0.5313 1.2063 

1.6724 

0.9011 

0.9200 

0.4000 0.9680 

1.4880 

0.8480 

1.4000 

5 1 0 0 0 




0 0 1 1 0 




0 0 0 0 1 . 




0.0000 167.20 

593.70 

156.50 

859.40 

-146.8 0.0000 

917.70 

104.40 

5695.0 

-78.31 -73.87 ( 

1.0000 

66.560 

634.80 

-114.8 -115.7 - 

472.60 

0.0000 

1300.0 

372.80 203.70 

-171.8 

342.40 

0.0000 

Five UNIFAC 

groups identified as : 

- ACH, - ACCIh, ■ 

- CH 3 CO, - CH 3 and - H 2 0 

Data Source : 

Magnussen ' 

'^1 


Simulation Run s System 2 

Feed conditions for steady state run (Run A) and perturbations for dynamic state 
(Run B) are listed in Table D.2. For dynamic simulation initial feed conditions has 
been taken corresponding to Run A at time t=0 seconds. 


Table D.2; Simulation Runs System 2 


Run 

Remark 

Feed condition/Perturbation 


Steady State 

Fig D.l to D.6 

S=0.32 moles/s, F=0.65 moles/s, S/F=2.3 w/w 

Aq. feed cone. (0.90, 0.10, 0.00) 

B. 

Dynamic 

Fig D.7 to D.ll 

(1) t=0 sec- feed condition as in Run A. 

(2) t=30 sec- F increased from 0.65 to 0.75 moles/s 

(3) t=400 sec- S increased from 0.32 to 0.42 moles/s 

(4) t,=80Q sec- Aqueous feed concentration changes from 
(0.90,0.10,0.0) to (0.8, 0.2, 0.0). 

(5) t— 1600 sec- Unperturbed all three perturbation 
(above) given so far. 

Allow simulation for 1 hour. 

''a IT ""4 i • t i r\r\\ 


S=Solvent phase (Toluene) feed rate (Average molecular weight = 92) 
F~ Aqueous phase feed rate (Average molecular weight = 20) 
S/F= Solvent to feed ratio 





MOLE FRACTION MOLE FRACTION 









Figure D.4: Forward flow (continuous aqueous phase) (Run A) 
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Stage l | 

Stage 2 
Stage 3 


TIME(sec) 


Figure D 

(Run B) 
0,2 


7: Transient profile for Acetone concentration in dispersed (Toluene) phase 


TJME(sec) 


Figure D.8: Transient profile for Acetone concentration in continuous (Aqueous) 
phase (Run B) 
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Figure I). 11: 'IVansient profile for forward flow of dispersed (Toluene) phase (Run B) 




Appendix E 

Elements of Jacobian 


This section gives elements of Jacobian shown in Figure 3.1. The Jacobian is highly 
sparse in nature so only non-zero elements are written in this section. Arrangements 
of the elements for any matrix block (A,B or C) has been shown in Figure 3.2. 
Each element represents the partial derivative of model equation (shown as rows in 
Figure 3.2) with respect to independent variable (shown as column in Figure 3.2). 
Partial derivative which have not been given in this section are zero. 
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